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A main feature of high-energy scattering in QCD is saturation in the number density of gluons. 
This phenomenon is described by non-linear evolution equations, JIMWLK and BK, which have 
been derived at leading logarithmic accuracy. In this paper we generalize this framework to 
include running coupling corrections to the evolution kernel. We develop a dispersive represen- 
tation of the dressed gluon propagator in the background of Weiszacker Williams fields and use 
it to compute C(/3q _1 q;™) corrections to the kernel to all orders in perturbation theory. The 
resummed kernels present infrared-renormalon ambiguities, which are indicative of the form and 
importance of non-perturbative power corrections. We investigate numerically the effect of the 
newly computed perturbative corrections as well as the power corrections on the evolution and 
find that at present energies they are both significant. 



1 Introduction 



Modern hadron collider experiments such as HERA, RHIC and especially the forthcoming LHC 
operate at high enough energies to observe new phenomenon associated with high gluon density. 
The principal characteristics of high-energy QCD scattering are the following: firstly, owing to 
Lorentz contraction, the configurations probed appear to be frozen in time compared to the natural 
time scales of the interaction. Secondly, the number density of soft gluons gets saturated at densities 
of the order of l/g 2 - These features, which are usually referred to as the Color Glass Condensate 
(CGC), are a consequence of the non-Abelian nature of the interaction and of the fact that gluons 
are massless. These features are therefore unique to QCD. 



1 



At sufficiently high energy the dominant interaction between the projectile and the target can 
be described by ensembles of boost-enhanced field configurations, the "frozen" modes mentioned 
above. This has been extensively explored in the context of the McLerran-Venugopalan model [1-9] . 
This description is tailored for asymmetric situations in which the field of, say, the target can be 
argued to be much stronger than that of the projectile. High-energy scattering of a virtual photon 
on a large nucleus is the prototype example of this situation. Nonetheless, at sufficiently high energy 
this generic picture is applicable to nucleus-nucleus scattering as well. 

Energy dependence can be incorporated into this picture by taking into account fluctuations that 
acquire properties of the previously frozen modes as one increases the collision energy. The relevant 
contributions are characterized by large logarithms ln(s) in the total invariant energy s in the 
collision. At low gluon densities, or weak fields, the resummation of high-energy logarithms has 
been formulated long ago as a linear evolution equation for the gluon distribution function, the 
BFKL equation [10-14]. However, at high densities the resummation of these logarithms leads 
instead to non-linear evolution equations for gauge field correlators. These can be formulated as 
a functional evolution equation known as the JIMWLK equation [6, 15-22], or equivalently, as 
an infinite coupled hierarchy of evolution equations for field correlators known as the Balitsky 
hierarchy [23-25]. A truncation of this hierarchy that retains most of its essential properties is 
known as the BK equation [23-27]. This equation describes high-gluon-density dynamics in terms 
of dipole degrees of freedom. In Refs. [26, 27] the BK equation has been derived from Mueller's 
dipole model [28-31], using nuclear enhancement as a tool to trace the dominant field configurations. 
This extends the ideas of the dipole model beyond the mere onset of saturation effects [32] . 

The most prominent feature of the solutions of these non-linear evolution equations is the emer- 
gence of an energy-dependent transverse correlation length i? s , or saturation scale Q s ~ ^/R s , 
which, asymptotically, encodes all the energy dependence of the cross section. The saturation scale 
characterizes the transverse momentum scales of radiated gluons that contribute to the evolution 
at any given energy. Modes much softer than Q s decouple: the number densities of soft gluons are 
saturated, so they remain constant as the energy increases. Independently of the initial condition, 
at sufficiently high energies the saturation scale increases rapidly with the energy. Therefore, Q s 
can be considered a hard scale: Q s 3> A. 

The possibility to describe saturation by perturbative evolution equations is a highly non-trivial 
result, since the gauge field involved is necessarily strong. The evolution equation is derived per- 
turbatively by expanding in small fluctuations on a strong background field. This is justified a 
posteriori: having found that soft modes do not contribute to the evolution, the equation is per- 
turbatively consistent. This kind of infrared stability is a direct consequence of gluon saturation. 
It is therefore not shared by the linear BFKL equation, which is instead afflicted by diffusion into 
the infrared. As was beautifully illustrated in Ref. [33] it is the non-linearity of the JIMWLK and 
BK equations that makes them infrared stable. The presence of the nonlinearities and hence Q s 
will also modify the role and influence of power corrections compared to the BFKL case discussed 
in Ref. [34]. 

Despite these strengths, JIMWLK and BK evolution suffer from a serious shortcoming: they are 
derived only at leading logarithmic accuracy, i.e. at fixed coupling. To partially compensate for 
this, all recent studies of the evolution have included running-coupling effects in some more or less 
ad hoc manner. There are several reasons why running-coupling effects are essential: 

• Running-coupling effects are known [35-43] to provide a large part of the next-to-leading- 
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order (NLO) corrections to the evolution in the low density limit, where the description 
matches onto the BFKL equation. 

• On the purely phenomenological side, for example in fits to HERA data [44], running coupling 
(or more precisely a dependence of the coupling on a scale involved in a single emission step, 
see below) is essential to slow down evolution by reducing gluon emission from small objects. 

• Conceptually, it is understood that the evolution is dominated by scales of the order of Q s . 
The non-linearity of the equation ensures that dipoles much larger than 1/Q S are inherently 
suppressed through the evolution (the saturation mechanism). However, with strictly fixed 
coupling, dipoles much smaller than 1/Q S still contribute to the evolution. As soon as the 
coupling depends on the size of the emitting dipolc such contributions are also suppressed 
through the evolution [45]. 

Despite both the practical and conceptual importance of running-coupling effects in the nonlinear 
evolution equations, there has been no derivation of how they enter. All simulations done so far 
involved ad hoc prescriptions for the scale of the coupling, based on nothing more than educated 
guesswork. 

In this paper we approach the problem on the more fundamental level. We show that the JIMWLK 
and BK equations can indeed be derived beyond the fixed coupling level. We find that the equations 
take a similar form to the fixed coupling case, while their kernel changes in a rather drastic way. 
For example, it does not naturally appear as a single scale-dependent coupling times the LO scale- 
invariant kernel. We explicitly compute running coupling C(/3Q _1 a") corrections to the kernel 
to all orders and resum them by means of Borel summation. We find that the resummed kernel 
present infrarcd-rcnormalon ambiguities. These arc indicative of the form and importance of non- 
pcrturbative power corrections. 

In order to perform this calculation we develop a dispersive representation of the dressed gluon 
propagator in the background of Weiszacker- Williams fields. This is a generalization of the well- 
known dispersive representation of the free dressed gluon propagator, a technique that has been used 
to compute running coupling corrections and estimate power corrections in a variety of applications, 
see e.g. [46-55]. 

As in the BFKL case, the non-linear evolution equations are expected to receive additional sub- 
leading corrections, which are not related to the running of the coupling. In this paper wc will not 
attempt to include such corrections. Some steps in this direction on the level of the BK equation 
have been taken in Ref. [56], or with an entirely different focus in Rcf. [57], and can be combined 
with our treatment where desired. 

The structure of the paper is as follows: in Sec. we give a short introduction to the physics 
described by the JIMWLK equation, in order to establish the key ideas and the notation that will 
be used in the rest of the paper. For more detailed background we refer the reader to the original 
literature [6,15-27] or review articles [58-61]. In Sec. 12.11 wc briefly review the JIMWLK and BK 
equations and recall those details of their derivation that are needed in what follows. Sec. 12.21 
is devoted to a discussion of running-coupling effects in JIMWLK and BK evolution, contrasting 
what has been done previously with what we want to achieve in this paper. This will be important 
also to clarify the terminology used in the remainder of the paper. Sec. |3| extends the derivation of 
the JIMWLK equation to the running-coupling case. It is divided into four subsections: Sec. 13.11 
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collects the tools of the conventional dispersive technique for the calculation of running-coupling 
corrections in the free field case. Sec. 13.21 generalizes these tools for use in the presence of the 
Weizsacker- Williams background as needed in the derivation of the JIMWLK equation. Sec. 13.31 
presents a re-derivation of the JIMWLK equation with a running coupling. Next, in Sec. 13.41 we 
formulate the newly computed corrections to the kernel as an all-order Borel sum. In Sec. 0] we 
discuss the convergence of perturbation theory. In Sec. [S] we use the renormalon singularities to 
determine the parametric form and the typical magnitude of power corrections affecting the kernel. 
Finally, in Sec. [5] we investigate numerically the effect of running coupling and power corrections 
on the BK evolution as a function of the saturation scale. In Sec. we summarize our conclusions. 

2 The physics of the JIMWLK equation 

The key points can be most easily understood in the context of deep inelastic scattering (DIS) 
of leptons on protons or nuclei, where q and p are the momenta of the virtual photon and the 
target respectively. Here, two kinematic variables play a role: (1) the deeply-spacelike momentum 
q 2 = —Q 2 < carried by the exchanged photon. This scale defines the transverse resolution of 
the probe and thereby the apparent size of the quarks and gluons encountered; and (2) Bjorken 
x := Q 2 /(2p ■ q), which is inversely proportional to the total energy s = (p + q) 2 in the collision: 
x w Q 2 /s. At high energy, the rapidity Y is directly related to Bjorken x via Y = ln(l/x). The 
rapidity is the natural evolution variable since ln(l/a;) « ln(s/Q 2 ) reflects the large hierarchy of 
scales in the high-energy limit, which appears with increasing powers in perturbation theory. 

At large Q 2 with fixed x there are well-established methods to treat such a system based on the 
Operator Product Expansion (OPE), a short-distance expansion in powers of l/Q 2 - Since Q 2 also 
controls the apparent size of the particles encountered, the OPE can be viewed as a small density 
expansion, despite the fact that particle numbers, driven by large logarithms in Q 2 , increase in par- 
allel with increasing resolution. As a consequence the description at large Q 2 can be based entirely 
on single-particle properties such as quark and gluon distribution functions: particle correlations 
are not important. This restriction is a key ingredient of the derivation of the Dokshitzer-Gribov- 
Lipatov-Altarclli-Parisi (DGLAP) equations that describe the increase of particle numbers with 
Q 2 in this domain. 

Going to small x at fixed Q 2 , no matter how large, one ends up in an entirely different domain, 
that of high densities, even if one starts out in a dilute situation. As the energy increases BFKL 
evolution keeps generating new particles (mostly gluons) which are all of effective size of ©(l/Q 2 ), 
and so the density keeps increasing. Eventually the density reaches a level where particle correla- 
tions become essential [28-32] and a description in terms of distribution functions alone becomes 
untenable. Since the BFKL description is based on gluon distributions, this is also the point where 
this evolution equation ceases to be adequate. Appropriate degrees of freedom and more general 
evolution equations are needed to describe the system beyond this point. The most general of these 
existing to date arc the JIMWLK equation, or, the completely equivalent Balitsky hierarchies, with 
their factorized truncation, the BK equation. 

JIMWLK and BK equations are formulated in terms of path-ordered exponentials, as defined in 
Eq. @ below, with paths collinear to the projectile direction, which can be interpreted as quark 
and gluon constituents of the projectile. The path-ordered exponentials encode the fact that, owing 
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to the high energy in the collision, these constituents penetrate the target without being deflected 
from their straight-line trajectories. The j*A cross section then reads 

W^,Q 2 )=Im™<cQ>™ = Jd 2 r\if(r 2 Q 2 ) J d 2 b / ^ ~"* U l) \ (y) (1) 

where r = x — y corresponds to the transverse size of a given qq dipole and b = (x + y)/2 
to the impact parameter of this dipole relative to the target, q 2 := — Q 2 is the large spacclikc 
momentum carried by the virtual photon. The square of the qq component of the photon wave 
function \ip\ 2 {r 2 Q 2 ) describes the probability to find a qq pair of size r inside the virtual photon 
and can be calculated in QED. It consists of a known combination of Bessel functions together with 
an integral over longitudinal momentum fractions already absorbed in the notation. The remaining 
factor, the expectation value of [/-operators is usually called the dipole cross section <7di P oie(X, r) of 
the target in question. All the properties of this interaction — details of the target wave function, 
gluon exchange between the target and projectile etc. — are encoded in this expectation value. 
The leading logarithmic corrections at small x appear in powers of a s ln(l/x). These corrections 
are resummcd by the JIMWLK equation. 

The dipole operator N xy := tr(l — U x ll y ) / N c itself is naturally bounded between zero and one. 
Typically, gluon densities grow towards large r, such that the expectation value of N xy interpo- 
lates between for infinitcsimally small dipoles, and 1 for very large ones. This encodes the idea 
of color transparency at short distances and saturation at large distances where gluon densities 
grow up to 0(1/ g 2 ). The length scale that characterizes the transition between the two domains 
can be interpreted as the correlation length R s of {/-operators, or equivalently gluon fields. The 
corresponding momentum scale Q s oc l/R s is usually called the saturation scale. Clearly, as more 
gluons are generated in JIMWLK evolution towards small x, the correlation length gets small and 
Q s increases. 

One key feature that emerges for this evolution is that details about the initial conditions are erased 
quickly and a universal scaling form of correlators such as the dipole cross section is reached. From 
then on, all xorY dependence is carried by the saturation scale Q S (Y). Such behavior has been seen 
in HERA data (geometric scaling) [62-64] and has important consequences for RHIC (disappearance 
of Cronin enhancement from mid to forward rapidities) [65-69] and the LHC experiments where 
the energies are higher. 

It had been noted early on in the context of the BK equation that a treatment at the strictly 
leading-logarithmic level is insufficient: running-coupling effects have a strong influence on the 
speed of evolution; quantities like the evolution rate X(Y) :— dy In Q 2 {Y) are reduced by more than 
50% if running-coupling effects arc introduced (in some heuristic way). Despite the explicit scale 
breaking introduced by the appearance of Aqcd in the running coupling, scaling of the dipole cross 
section with Q s is retained to very good accuracy. In Ref. [45] it was emphasized that running- 
coupling effects are also conceptually important: only with running-coupling effects included does 
the phase-space region active in the evolution center around the physical scale Q s . At strictly 
leading-logarithmic level — in the conformal limit for the evoluton kernel — the evolution involves 
short-distance 1 contributions from more than 7 orders of magnitude away from the physical scale. 

J As explained above, the infrared is not a problem due to the presence of the correlation length R a , which acts 
an an effective infrared cutoff. 
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Unlike the other evolution equations in QCD, such as DGLAP and BFKL, the JIMWLK and BK 
equations have been derived only at leading logarithmic accuracy. Only partial calculations of 
two-loop corrections to the BK equation are available [56] but they do not include any attempt to 
determine the running of the coupling. For our purposes, the existing results in the low-density 
limit for BFKL evolution are of limited use: they offer no hint as for how to extend or extrapolate 
the calculation into the high-density domain where non-linearities appear. A direct calculation in 
the context of the JIMWLK and BK equations is therefore necessary. 

As announced in the introduction, in this paper we compute running-coupling corrections to the 
JIMWLK and BK equations. In Sec. EH we briefly review the fixed-coupling derivation of the 
JIMWLK equation, by considering small fluctuations in a strong Weizsacker- Williams field that is 
encoded in the eikonal factors of Eq. . This derivation will then be generalized to the running- 
coupling case in Sec. |3] using a dispersive representation of the dressed gluon propagator in such a 
background field. This will not only enable us to calculate running-coupling corrections, but also 
to explore non-perturbative effects in the evolution. 

2.1 Evolution equations 

The JIMWLK equation and the Balitsky hierarchies 

The JIMWLK equation is a functional Fokker-Planck equation for the statistical weight Zy\U] 
defining the y-dependent averaging procedure (. . .)(Y) that determines the expectation value of 
operators 0[U] made of an arbitrary number of path-ordered exponentials U, where 



Such an average was already encountered above in Eq. (JIJ in the case of the dipole operator 
tv{l — U x U^)/N c . To understand why these averages determine virtually all cross sections at 
small x, recall the origin of the path-ordered exponentials J2J): they encode the interaction of 
fast moving quarks and gluons in the projectile wave function with the target field. It is the high 
energy of the collision that allows the description of this interaction in terms of the eikonal factors 
U, that at leading order follow perfectly straight, lightlike worldlines. We have chosen a frame 
where these trajectories extend along the minus light-cone direction at x + = 0. Each particle is 
then characterized by the remaining coordinates, namely its transverse location x. The leading con- 
tribution comes from interaction with the non-Abelian Weizsacker- Williams field of the target, A + . 
It takes the form 



where 6 + , or more specifically (3(x), is the single leading degree-of-frcedom, a strong field, while 
5 A is a small fluctuation in which we will expand. The S(x~) reflects the lack of resolution in 
longitudinal direction: no internal details of the field of the target are probed. The independence 
on x + reflects the fact that the target wave function is frozen during the interaction, an extreme 
time dilation. At fixed rapidity all dominant contributions are determined by the background field 
b + = 6(x~)f3(x). Moreover, they only enter in a very specific form, via the Wilson lines U. 




A+{x) = b + ( x -,x) + SA{x), 



b+( x -,x) = S(x-)(3(x) 



(3) 
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For some generic operator made of these Wilson-line fields, 0[?7], the average of Eq. JTJ will be 
written as 

(0[U]) Y := J b[U]0[U]Z Y [U] (4) 

where D[U] is a functional Haar-measure and Zy [U] contains the detailed physics beyond the cikonal 
approximation already incorporated by selecting U as the relevant degrees-of-freedom: Zy[U] is the 
statistical weight for all possible field configurations. 

Most of our knowledge of Zy[U] is perturbativc: the JIMWLK equation 2 determines the Y depen- 
dence of this average: 

dyZy[U] = -H[U]Zy[U], (5) 

where the JIMWLK Hamiltonian TC[U] is given by 

H[U] = -^tp-JC^y [Uf (iV a JV b y + iVS*V») + (iV a JK + iV£*V2) } , (6) 
with the LO kernel: 



r — v ~ ; v ~ y > — ' i ' 2 f7 s 

^ XZ V _ t~ _ _ „\2 ~~ ^2^2' v> 



(g - z) ■ (z - y) ri -r 2 

(cc - z) 2 (z - y) 2 r 1 2 r 2 ' 2 

where x, y and z are transverse coordinates. Here we also introduced a shorthand notation for the 
vectors connecting the points in the transverse plain: 



r = x-y . n = x-z , r 2 = y- z , (8) 

and their lengths: r± := |t"i|, etc. The notation in © assumes an integration convention over 
repeated coordinates appearing as an index in /C xzy and in the vector field operators V£. The 
Hamiltonian H[U] is second order in left- and right-invariant vector fields V° and V°, which are 
Lie derivatives: they act on the Wilson-line variables U x according to 3 

iV a x U y := -U x t a 5W ; iV a JJ y := t a U x 5^ . (9) 

The terms in the square brackets on the r.h.s of Eq. © arc grouped according to their origin in 
real-emission and virtual corrections, respectively: real-emission contributions involve an additional 
Wilson line U^ b at transverse location z. 

The full derivation of Eq. JSJ has been presented exhaustively in Rcfs. [20-23]. We nevertheless 
need to recall here how the JIMWLK Hamiltonian relates to Feynman diagrams, in order to prepare 
its re-derivation with running coupling corrections. Rapidity dependence of Zy[U] at LO is driven 
by the lowest-order fluctuations 5A around the background 5(x~)(3(x) of Eq. 0. In this sense, 
the LO JIMWLK Hamiltonian in Eq. JfJJ) is constructed such that it adds the LO "exchange" and 
"self-energy" corrections to, say, an interacting qq pair, represented by its Wilson-line bilinear 

2 For a first derivation of this equation in its most compact form see Ref. [20]. The version presented here is based 
on Ref. [60]. 

3 See Ref. [20] for more details. 
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u x ®W y . 



ln(l/a) H[U] U x ®Ul = 




The diagrams shown in (|10|l are Feynman diagrams where the gluon propagator of the fluctuations 
(5 ASA) is taken in the background of the strong target field S(x~)(3(x). 

The correspondence to real and virtual diagrams becomes visible upon resolving the Feynman 
diagrams into x~ ordered diagrams of light-cone perturbation theory For instance, 




where the diagrams on the r.h.s. should be interpreted as diagrams of light-cone perturbation 
theory. Light-cone time x~ runs from bottom right to top left, the two collinear 4 Wilson lines in 
this direction represent the dipole (the projectile) and the target is shown as a perpendicular line at 
xr = 0, from bottom left to top right. To be precise, the last two diagrams in are a shorthand 
notation for a sum of two different x~ orderings each: 




(12) 

The factors of U, U> and U ab representing the interactions of a projectile quark, antiquark and 
gluon, respectively, are indicated as large dots where these Wilson lines cross the target line. In 
a derivation of JIMWLK based on projectile wave functions 5 , the contribution of the first two 
diagrams of (| 1 1 1> are associated with a situation where the interacting gluon reaches the final state 
and in this sense they correspond to real-emission diagrams. Correspondingly, as indicated by 
the third dot, they contain an additional Wilson line for this produced, eikonally interacting gluon, 
which appears in ® as Uf. The remaining diagrams represent purely virtual contributions in which 
the number of Wilson lines does not change. Similarly, both real emission and virtual corrections 
are present in the self-energy-like diagrams in (|10f) . 

4 Notc that these two lines are actually separated only in the transverse direction. 

5 See Refs. [28-31] for a derivation of JIMWLK based on projectile wave functions in the context of the dipole 
model, Refs. [26,27] for the BK case, or Ref. [70] for a re-derivation of JIMWLK via amplitudes. 
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Note that in the absence of the target field the eikonal factors become trivial: U — ► 1, and then 
there is strictly no evolution. In this limit each individual diagram on the r.h.s. of (|10|) vanishes 
identically owing to exact real-virtual cancellation: in the first two diagrams, corresponding 
to real gluon emission, cancel against the last two diagrams, which represent virtual corrections. 
Analogous cancellations occur in this limit in the light-cone perturbation theory decomposition of 
the self energy like diagrams in JTUI) . 

As a functional equation, Eq. J^J) is equivalent to an infinite set of equations for n-point correlators 
of U and W fields in any representation of SU(N C ), called the Balitsky hierarchies [23]. In order to 
obtain the evolution equation for a correlator of a given composite operator 0[U] made of [/-fields, 
one first takes the Y derivative of Q and then uses JSJl to replace the dyZy[U] by —H[U] Zy[U], 
obtaining: 

d Y (0[U}) ( Y ) = -[ D[U] 0[U] H[U] Z Y [U], (13) 

where the Hamiltonian still acts on Zy[U]. Using the self-adjoint nature of Tt[U] with respect to 
the Haar measure (as ensured by the Lie derivatives), one can then rewrite (|13|) as 

d Y {0[U})(Y) = -((H[U}0[U]))(Y). (14) 

Here the JIMWLK Hamiltonian acts on 0[U] . Observing that H[U] explicitly contains a U -operator, 
and that the number of [/-operators remains invariant when acted upon by the Lie derivatives V° 
(see Eq. one understands that the evolution equation for (0[U\) must involve operators with 
more U fields than the original operator 0[U}. Thus, the nonlinear nature oiH[U] implies that the 
r.h.s. of the equation depends on a new type of correlator of U fields. To determine the rapidity 
dependence of (0[U])(Y) one will therefore need also the evolution equation of this new correlator, 
which in turn will couple to yet higher composite operators containing more U fields. Continuing 
the process one ends up with an infinite hierarchy of equations, defined by the operator 0[U] used 
to start the process. The derivation of the BK equation shown below provides a simple example 
for such a hierarchy. 



Truncation and the BK equation 

As the simplest correlator with immediate phenomenological relevance we consider the two-point 
function of the dipole operator 6 : 

Ny, xy :={N xy ){Y) N xy := tT{1 ~^ Uy) . (15) 

Using Eq. (|14|) in the case of N xy with the explicit expressions corresponding to the diagrams in 
Eq. (|l(Jfl . one immediately obtains: 

dy(N xy )(Y) = Jd 2 z(2JC xzy - K xzx - JC yzy ){N xz + N zy - N xy - N xz N zy ){Y). (16) 

The linear combination of JIMWLK kernels K, that appear in this expression are in one-to-one 
correspondence with the three diagrams in . They combine into the very compact form of the 

"Below we will often use a shorthand notation where the Y dependence is not explicitly indicated. 
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BK kernel IC xzy : 



(x-y) 2 
(x - z) 2 (z - y) 2 



K^xzy - ^^-xzy ^xzx K-"yzy ~ V2 ' (^^) 



Clearly the r.h.s. of Eq. ijltjfl depends on a 3-point function containing operators with up to four 
f/W factors, that in general does not factorize into a product of two 2-point correlators: 

(N xz N zy )(Y) = (N XZ ){Y) x (N ZV ){Y) + corrections . (18) 

To completely specify the evolution of (N xy )(Y) one therefore need to know (N xz N zy ){Y). The 
latter, in its evolution equation, will couple to yet higher n-point functions and thus one is faced 
with an infinite hierarchy of evolution equations, as anticipated above. The entire hierarchy (as 
well as others, corresponding to other composite operators 0[U}) is encoded in the single functional 
equation JSJ). If one drops the corrections in l|ltj[) and factorizes the correlators in the spirit of a large- 
N c approximation, the hierarchy is truncated and reduces to a single equation. This truncation can 
be interpreted as an independent scattering approximation and it leads to the Balitsky-Kovchegov 
(BK) equation [24,27]: 

d Y N Yxy = ° s ^J Nc I d 2 z jC xzy [n Yxz + N Y , zy - N Yxy - N Yxz N Yzy ) . (19) 



Generic features and infrared stability 



Despite the complex nature of the evolution equation, it is possible to gain insight into some generic 
features. It can be proven [20] that the evolution equation possesses an attractive fixed point at 
Y — > oo at which the system has vanishing correlation length. For the evolution of a physical 
correlator, such as N Y ^ xy , this implies a generic trend as shown in Fig. Q](a): with increasing Y 
saturation (N YtXy — > 1) is reached at ever shorter distances. This leads to a further important 
property of the evolution equation, namely its infrared stability: the evolution is not affected 
by long-wavelength fluctuations beyond the characteristic correlation length (the inverse of the 
saturation scale Q s ) where N Y xy = 1. In this way the saturation scale acts as an effective infrared 
cutoff. Fig. n (b) demonstrates that the modes that contribute to the evolution have momenta of 
order of Q S (Y). With increasing Y, the active region moves towards the ultraviolet. 



2.2 What we mean by running coupling 

Similarly to other evolution equations in QCD [71], the LO JIMWLK kernel is conformally invari- 
ant 7 . As usual, breaking of this symmetry is expected to appear through radiative corrections at 
0(/3oa 2 ) where (3q is the leading coefficient of the f3 function, 

da ^I n - -0o Mm 2 )/-) 2 - (h Mm 2 )/-) 3 + • • • ; A) = %c A l -N f . (20) 

7 Note that while the LO kernel is conformally invariant, the solution of the equation is not: it is characterized 
by a correlation length. This scale originates in the initial condition, and it is preserved owing to the non-linearity 
of the equation. 
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Figure 1: Generic evolution trend for a single-scale dipole correlator, (a) shows N(Y,r) as a function of 
the dipole size r = \x — y\ for several values of the saturation scale Q S (Y). Q s increases with Y ; saturation 
then sets in at smaller distances, (b) shows dyN and thus the activity in a given evolution step as a function 
of the same variables. With increasing Q S (Y) contributions are centered at ever shorter distances. 

At higher orders one expects corrections C(/3Q _1 a"), as well as ones associated with subleading 
coefficients of the (3 function, e.g. 0((3\(3q ~ 3 a"). In Sec. |3]we compute these running-coupling 
corrections using a dispersive representation of the dressed gluon propagator in a background field. 
We will show that while the general structure of the evolution equation, namely Eq. ©, holds 
as at LO, the kernel itself changes drastically with respect to Eq. Q — see Eq. (|62fl or l|77[l 
below. A similar generalization holds in the BK case, see Eqs. (|91|) and i|92[l . JIMWLK and BK 
evolution with running-coupling is therefore qualitatively different from the fixed-coupling case in 
that at each step in the evolution the coupling depends on the details of the evolving configuration; 
moreover, different final states, that are characterized by different "daughter dipole" sizes, are 
weighted differently in the r.h.s. of the evolution equation. 

Although all previous derivations of JIMWLK and BK evolution equations were restricted to the 
leading logarithmic approximation, it has been clear for quite a while that running-coupling correc- 
tions will be necessary to get quantitative results 8 . By making physically-motivated scale choices 
for the renormalization scale of the coupling, several authors found that the evolution rate X(Y) 
reduces by a factor of two or more compared to the fixed-coupling case. 

It is important to make a clear distinction between the actual higher-order contributions to the 
kernel, which we compute in the following, and ad hoc choices of scale for the coupling, which 
have been often referred to as "running coupling" . Such prescriptions have been assumed in all 
numerical simulation of the BK equation [33,45,72-77]. Similar assumptions were used in the 
analytical estimates of Ref. [78], which agree well with numerical simulations [45]. To put the 
results of the present paper in context of previous work, we find it useful to further distinguish 

8 In Ref. [45] it has been further emphasized that running-coupling corrections, where the scale of the coupling 
depends on the scales involved in a single evolution step, are needed to reduce active phase space in the ultraviolet 
from about 6 orders of magnitude to about one. Despite this dramatic reduction of phase space, other qualitative 
features of the solution of JIMWLK or BK equations were found to be the same: for any set of initial conditions that 
interpolate between color transparency at short distance and saturation at large distance (as shown schematically in 
Pig. the system approaches an asymptotic line where the dipole correlator reaches a near— scaling form. 
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between different categories of scale choices that were made in the literature: 

1. Fixed or essentially fixed coupling: a s (^ 2 ) is treated as a constant or as a function of Y. 
Within this class, simulation results, starting from the same initial condition, can be related 
by re-scaling the evolution variable. For example, the scenario where the coupling depends on 
Q S (Y) can be related to the one where the scale of the coupling is fixed as some Qo through 
the change of variables (see also [79]): 

_ a s (QgfY)) 

«. (Qo) 

Obviously, any ansatz of this sort would fail to capture the essential physics of running 
coupling, and would inherit the ultraviolate phase-space problem discussed in Ref. [45]. 

A non-trivial change in the shape of the solutions with respect to the fixed-coupling case occurs 
if the scale of the coupling is determined by the size of the "dipoles" involved in the evolution. 
In such cases the ultraviolate phase-space problem is generally removed, since emission from very 
small objects is suppressed by small coupling constants. The remaining two cases fall into this 
category. 

2. "Parent dipole" running: where the scale of the coupling on r.h.s. of the BK equation l|19|l 
is assumed to depend on the initial dipole size r = \x — y\, namely 

a s (p 2 )JC xzy -> a s (c 2 /r 2 )/C xzy , (21) 

where c is a dimensionless "scale factor" of order one 9 . This has been the most common 
ansatz in numerical simulations of the BK equation. 

While this appears to be a natural ansatz in the context of the BK equation, it can not 
be easily reconciled with the JIMWLK equation. To see this, recall that © is a functional 
equation. Therefore, the Hamiltonian in Eq. 10 cannot depend on any scales characterizing 
a particular operator 0[U]; all three transverse coordinates appearing in (JSJ are integrated 
over internally. When deriving the BK equation from JIMWLK as in (|16|l one obtains the 
BK kernel as a specific combination of the original JIMWLK kernel: 

OtsfcxZy * ^ tt s (fl xzy )lC XZ y ~ Oi S (^ xzx ) JC XZX — U S (^y Zy ) K.y Z y, (22) 

where each term corresponds to one of the diagrams in HUfl . Whatever one assumes of the 
functional dependence of n xzy on the coordinates, since the self-energy-like diagrams are 
independent of the parent dipole size \x — y\, Eq. (1221) does not lead to Eq. (|21|l . 

Even considering the BK case on its own, "parent-dipole" running may appear artificial: there is 
no reason why the sizes of the dipoles produced in the same step in the evolution, which can of 
course be quite different from the parent size, should not play an important role [77]. This leads us 
to the final category: 

9 It was often taken as c fa 4, motivated by direct Fourier transform in the double logarithmic limit. 
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3. Final— state— dependent evolution: where the scale of the coupling on r.h.s. of the BK 
<|19|) or JIMWLK (JSJ equations is assumed to depend on all three distance scales present in 
a single evolution step, namely the "parent dipole" r = \x — y\ and the two newly produced 
"daughter dipoles" , n = \x — z\ and r2 = \y — z\. Obviously, in this case the coupling affects 
the weight of different final states, depending on the transverse location of the emitted gluon 
(the coordinate z) that is integrated over. 

Several such models were proposed, e.g. \/a s (c 2 /rf) y / a s (c 2 /r§) or a s (max {c 2 /r 2 , c 2 , 
and found to be in fair agreement with "parcnt-dipole" running, sec e.g. [77]. As already 
mentioned, no deep justification of any of these models has been provided. 

In the next section we compute running-coupling corrections to the JIMWLK equation. In Sec. 
shall use this result to derive the corrections for the BK case and study the consequences numerically. 
As we will see, these corrections do depend on the details of the final state at each evolution step 
and involve all three scales. Moreover, in neither of the two equations do these corrections naturally 
reduce to a single scale-dependent coupling times the LO kernel. 

3 Derivation of JIMWLK evolution with running coupling 

3.1 Running coupling, Borel transform and the dispersive approach 

We are interested in improving the leading-logarithmic result of the JIMWLK and BK equations by 
including running-coupling effects, namely corrections associated with the renormalization-group 
equation (|20[) . In general, running-coupling corrections are important in QCD, for two reasons: 

• They usually constitute a large part of the higher-order (notably NLO) corrections [48,49, 
52,80-84]. The main reasons for this are: (1) the average virtuality of gluons is usually 
different (typically much lower) than the principal hard scale, which is often used as the 
default renormalization point; and (2) flo is sizeable. These effects are especially important 
when the hard scale is low, since then the coupling is large and its evolution 12(J|) is fast. 

• They dictate the large-order asymptotic behavior of the perturbative series, which is domi- 
nated by factorially increasing contributions, O (n!/3p _1 (a s /7r)") , the renormalons [85,86]. In 
this way the resummation provides some insight into the non-perturbative side of the prob- 
lem: the ambiguity in summing the perturbative expansion, which is expected to cancel in 
the full theory, indicates the parametric dependence of non perturbative corrections on the 
hard scales (power-suppressed corrections) and provides some clue on the potential size of 
these corrections [87-90]; for a review see Refs. [91,92]. 

Both these aspects are relevant in non-linear evolution. As demonstrated in Sec. 12. 21 the multi-scale 
nature of the problem calls for a systematic study. 

In this section we wish to briefly recall some basic ideas and techniques for resummation of running- 
coupling corrections that will be generalized and applied to the JIMWLK case in what follows. To 
this end, consider some perturbatively calculable (infrared and collinear safe) quantity, R(Q 2 /A 2 ), 
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depending, for simplicity, on a single external scale Q 2 and having an expansion starting at order 
a s (/i 2 )/7r with the LO coefficient normalized 10 to cqq = 1: 



R(Q 2 /A 2 



Cll 



ClO 



a s (/u 2 



(23) 



where (3q is defined in i|20|) and the cio term is a conformal coefficient, not associated with the 
running coupling. We will work in the large-/?o limit, where the cto term is formally subleading. 
The well-known BLM prescription [80] absorbs the NLO contribution that is leading in /3 into the 
LO by a scale choice: fi 2 LM = Q 2 exp{— en}. In this way also higher-order corrections often become 
smaller. This becomes intuitive upon looking at the momentum integral that is approximated by 
^(cLm)/ 1 ! where /U| LM acquires the interpretation of the average gluon virtuality. The all-order 
resummation of running-coupling effects in the single dressed gluon approximation, 



R(Q 2 /A 2 ) 





MQ 2 )\ 


1 large a 





1 



E 



00 



a s (Q 2 ) 



(24) 



can be viewed as a generalization of this procedure, where instead of a single optimal scale choice, the 
proper (observable-dependent) weight is given to any specific gluon virtuality. In the context of the 
dispersive approach presented below (see Eq. \6ty ). this weight function is called the "characteristic 
function" [50]. 

Technically, the resummation of all C(/3q _1 a") higher-order corrections becomes feasible in QCD, 
owing to its simple relation with the resummation of diagrams with an arbitrary number of fermion- 
loop insertions. Because Po is linear in Nf, one can simply compute the fermion-loop chain diagrams 
and then replace Nf by the non Abelian value of — 6/?o, according to ({20(1 . In more formal terms, one 
begins by considering the large- Nf limit with fixed Nfa s , the leading term in the flavor expansion. 
Clearly this limit itself is not physically interesting, it just provides a tool to identify running- 
coupling contributions in the approximation where the non- Abelian /3 function 120fl is one loop, the 
so-called large-/3o limit. 

The calculation of a gluon propagator, dressed by fermion-loop insertions is simplified by the fact 
that the fermion loop itself is transverse, 



iWfc 2 



k 2 Tl(k 2 ), 



and therefore, in any gauge such insertions affects only the propagating particle pole 1/k 2 . The 
all-order sum builds up a geometric series giving rise to a factor 1/(1 + H(k 2 )). The resummed 
propagator takes the form 11 : 



covariant gauges: 
strict axial gauge: 



1 



1 



k 2 1 
1 



n(fc 2 ) 
l 



k 2 i + n(fc 2 ) 



9V 
9»u 



1 kftki/ 



k 2 j s /t 2 k 2 

k^n v + n^k u 



k^kuTi 



k ■ n 



(k ■ n) 2 
[93] 



(25a) 
(25b) 



10 Using this normalization R(Q 2 / A 2 ) can also be interpreted as an "effective charge 1 

xl Note that the £ term in I25al is not affected as it does not describe a physical mode. £ is the width of a 
Gaussian approximation to a functional <5-function that is meant to implement the covariant gauge 9^71^ = 0. Only 
for vanishing width, £ = 0, do all gauge fields obey the gauge condition strictly. The axial propagator is written for 
vanishing width, i.e. the strict gauge condition. That is why the resummed terms multiply the whole structure. 
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with 

n<* 2 )L._, TO = ^ta(-* 2 e-!//), (26) 

where the renormalization of the fermion-loop n(fc 2 ) was done in the MS scheme 12 and where we 
already made the replacement: Nf — > — 6/3o- 

The next step would be to insert the dressed propagator into the relevant Feynman diagrams and 
perform the momentum integration, d A k. Since the sum over any number of II(fc 2 ) insertions has 
already been done, by performing the fc-integration one would hope to get directly the resummed 
physical quantity R(Q 2 /A 2 ) in Observing that the resummed propagator with has a 

Landau singularity, one realizes that this direct all-order calculation cannot be done. As we explain 
below, a rcgularization of the sum is required even if the coefficients c nn are all finite 13 . The simplest 
and most familiar way to see this is to represent the effective running coupling 14 that includes the 
dressing, as a Borel sum: 



aY(-k 2 -iO) _ a s (^) 1 _ _1_ 



T{u) (-£: 2 c-i/A 2 ) " , (27) 



7T 7T 1 + U(k 2 ) 0o J Q 

where for one-loop running coupling 15 , namely upon using l|2()|) . one simply gets 

a v s {-k 2 -iQ) 



one— loop 



Po ln(-fc 2 c-|/A 2 ) ' 



and therefore T(u) = 1. ft is now possible to proceed with the calculation of the Feynman diagrams 
where the only change is that the particle pole is modified into a cut: 

(28) 



-k 2 -io {-k 2 -ioy+ u ' 

This modification of the propagator is known as Borel or analytic rcgularization. 

The all-order resummation of a given quantity R in the large-/3o limit can therefore be done 
using l|27|) by first performing the momentum integration with the modified propagator i|28[l . which 
directly yields the Borel representation of the sum in the single dressed gluon approximation, 
namely: 

oo 

i?(Q 2 /A 2 )| largc/3o = 1 J duT(u) (Q 2 /A 2 y u B(u). (29) 

o 

Upon performing the momentum integration for a typical observable in QCD one would find that 
B(u) has singularities along the positive real axis, which is the integration axis. This obviously 

12 We denote the MS coupling a™ s {p?) by a s {p?). 

13 This is indeed the case if R is an observable, i.e. it requires no additional renormalization and has no infrared 
divergencies. 

14 It is convenient to absorb the factor e — 3 from the renormalization of the fermion loop in 1261 into the definition 
of the coupling. We follow this convention and define aY as the coupling in the V scheme, which is related to the 
MS scheme by Ay = es A 2 . 

lo Although we restrict the calculation to one— loop coupling where T(u) = 1, we will keep writing T(u) in any 
Borel representation: this allows for a straightforward generalization to two-loop running coupling following [94,95]. 
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means that Eq. (|29[) is ill-defined as it stands. The obstructing singularities are called infrared 
renormalons, and they merely reflect the fact that the series for R in 1|24|) and thus also in 1|23|) . is 
non summable. The problem only becomes manifest if one attempts to sum the series: any finite 
order expansion — here all the coefficients c nn in (|24J) — are well-defined and finite. They can be 
obtained by expanding B(u) under the integral, 

oo 

Bin) = ( 3 °) 

The Borcl function typically has a finite radius of convergence Uq, so c„„ grow as uI/uq at high 
orders 16 ; in the case of infrared renormalons in QCD, the singularity is at u = Uq > 0, so the 
coefficients all have the same sign 17 . In terms of the momentum integration these large perturbative 
coefficients arc associated with extremely small virtualities k 2 — > 0, namely soft modes whose 
dynamics is non-perturbative. Thus, the non-existence of the sum is a reflection of the fact that 
perturbation theory does not fully describe the dynamics, and the observable is sensitive to some 
extent to non-perturbative contributions. The ambiguity can therefore be resolved by a proper 
separation prescription between perturbative and non-perturbative corrections, either by an infrared 
curoff, or by other means, for example by modifying the integration contour of the Borel integral l|29[) 
in a particular way or by taking its principal value. The same prescription must then be applied to 
regularize the non-perturbative contribution. 

An additional, general property is that the singularities in B(u) occur at integer — and sometimes 
half integer — values of u. This corresponds to the fact that alternative definitions of the sum 
of the series (that arise e.g. by modifying the integration contour in (129(1 ) differ by integer - 
or half integer — powers of A 2 /Q 2 . These ambiguities must be cancelled by non-perturbative 
power corrections, and they can therefore serve as a perturbative probe of such effects. In cases 
where an operator product expansion applies, one can get a direct interpretation of the source 
of each ambiguity in terms of local operators of higher dimension (or twist) [87,89]. In such 
cases it is also possible to trace the cancellation of ambiguities [91,96,97]. In the absence of an 
operator product expansion, the rcnormalon technique often provides a unique window into the 
non-perturbative regime: by identifying the ambiguities in summing the perturbative series one 
learns about the parametric dependence of power corrections on the hard scales and about their 
potential size [46-55,87-91,96-99]. 

Computing directly the Borcl function B(u) using H28J) may sometimes present technical difficulties. 
One of the most effective techniques to deal with this problem is the dispersive technique, which 
has been used in a variety of applications, see e.g. [46-55,98,99]. In the next sections we shall 
generalize this technique to the case of Weizsacker- Williams background fields, in order to use it in 
the re-derivation of the JIMWLK equation. Let us therefore review here the basic idea, and collect 
the necessary formulae. 

The dispersive method, see e.g. [49,84,98], recasts the dressed gluon propagator, and thus the 
running coupling of (|27|l as a dispersive integral. As we shall see, this allows one to compute a 
generic quantity R(Q 2 /A 2 ) to all orders in the large-/?o limit, by simply replacing the massless 

16 This asymptotic behavior can be modified by additional factors of the form n 1 '^ , depending on the structure 
of the Borcl singularity. For simplicity wc assume here simple poles. 

17 In contrast, for ultraviolate renormalons in QCD uq < 0, so there is sign oscillation at high order. 



16 



gluon pole by a massive one, namely: 



-k 2 - iO 

instead of (l2*g|) . 

The dispersive representation takes the form 
aj(-k 2 -i0) a s ( M 2 ) 



i + n(fc 2 



*0' 



1 

Jo 



dm 



Pv(m 2 ) 
i 2 -k 2 - iO 



(31) 



(32) 



where is the discontinuity of the coupling on the time-like axis, 



Pv(m 2 ) 



Im{a^(- 



m 



i0)/7r} 



ip a s ( P 2 ) Im{n( ? 7i 2 )} 



7T |i+n(?7i 2 )| 2 

Alternatively, it is convenient to express Eq. I|32|) in terms of the "effective 18 time-like coupling 
a v s {-k 2 -i0) 



(33) 



= — I dm 

Po 



Pv(m 2 ) 



k 2 — iO /3o Jo m 



dm* 



dm 2 k 2 



(34) 



where A^ s {m 2 ) is defined by 



pv{m 2 ) =: m* 



dm 2 



A v cS {m 2 ). 



(35) 



The explicit expression for A^ s (m 2 ) in the case of a one-loop running-coupling (namely 13211 
with (HHJ) takes the form: 



A V oS {m 2 \ 



i 1 1 

= arctan 

lone— loop 2 7f 



1 m* 

7T A?, 



A two-loop expression for A^ s (m 2 ) can be found in Ref. [99]. 

The all-order rcsummation of R in the large-/?o limit can now be done using (|34fl by first performing 
the momentum integration with a massive gluon (i.e. with the replacement of Eq. I|31[ll leaving the 
dispersive integral over the mass undone. This leads to 19 



R{Q 2 /A 2 )\ 



1 



large f3 



dm* 



Po Jo 77)2 



^eV™ 2 ) 



m 2 — i T(m 2 IQ 2 ) 
dm* 



1 



dm* 



Po Jo m2 



pv(m 2 )[F(m 2 /Q 2 )-F(0)], 



(36) 



18 Lct us emphasize that we do not make any use of the dispersive representation to impose analyticity and regularize 
the Landau singularity, as done for example in Ref. [100]. We merely use the dispersive method to compute the 
Borel function, which serves as a generating function for the perturbative coefficients according to 1301 and as a tool 
to analyze infrared renormalon ambiguities. 

19 Note that we ignore pure power correction terms (A 2 /Q 2 ) which distinguish this integral from the Borel sum we 
eventually compute. These will be totally irrelevant here. The interested reader is referred to [49,51,52]. 
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where in the m 2 — ► limit J 7 (m 2 /Q 2 ) yields the LO coefficient, .F(O) = coo = 1. By definition, 
J r (m 2 /Q 2 ) depends on the specific quantity computed and it encodes perturbative as well as power- 
correction information (see e.g. [49,84]). T(m 2 /Q 2 ) is therefore called the "characteristic function" 
of the quantity or process under consideration. Upon expanding the effective coupling A^ s {m 2 ) 
in l|36|) in powers of ln(m 2 /A 2 ) and integrating term by term (log moments) one obtains the per- 
turbative coefficients c nn . On the other hand, expanding T(m 2 /Q 2 ) at small m 2 and extracting 
the non-analytic terms [84] one obtains information on renormalon ambiguities and thus on power 
corrections in A 2 /Q 2 . 

As explained above, the perturbative and power-correction information encoded in the characteris- 
tic function can be conveniently extracted from the Borcl representation (|29|) . The Borcl formulation 
also offers a natural definition of the sum, using the principal value (PV) prescription. The PV Borcl 
sum can be related to a Euclidian-cutoff regularization [52], and it has good analytic properties 
as a function of the hard scales [54, 101], e.g. it has no Landau singularities and it is real-valued. 
In order to relate the dispersive and the Borel formulations, jSSjl and iJSSJ, respectively, one first 
computes the Borel representation of the discontinuity of the coupling by using 127( 1 in 1|33|) : a 
straightforward calculation yields: 

/•OO 

p v (m 2 ) = - duT(u) S -^ei u (m 2 /A 2 y u . (37) 

Next, by substituting (|37|l in i|36|) and changing the order of integration one recovers the Borcl 
representation (|29|l with 

oo oo 

B{u) = _ e |«™ [dCC^WQ-m] = -e** — /dcc-^l. (38) 
TT J KU J dQ 

o o 

In the following wc will use this formula to compute the Borcl function B(u) from J-(C)- 

Finally, it should be emphasized that the resummed result for R(Q 2 /A 2 ), computed with a single 
dressed gluon, is exactly renormalization-scale invariant, despite using just a subset of all pertur- 
bative corrections. At power level, we shall define the sum using the principal value prescription 
of (J2I! and use the renormalon poles to analyze non-pcrturbative power corrections. 

3.2 The dressed propagator in the background field 

At this point, we have described the method used to derive the JIMWLK equation at fixed coupling 
fSec. l2.l|) and the generic tools used to promote a one-loop calculation to one with a dressed gluon 
fSec. 13. l|) . These need now to be combined. Since all the diagrams in Hl()|l involve a single gluon, 
the dispersive technique should be applicable. The crucial step is to construct a propagator that is 
both in the background of Weizsackcr- Williams field and is dressed by running -coupling corrections. 
This is the main goal of this section. 

Independently of the particular technique adopted, the calculation of running coupling corrections 
to the kernel involves a subtlety which is associated with the presence of new production channels, 
namely the production of a well separated qq pair or a pair of gluons, starting at NLO. These 
processes involve two additional Wilson lines — not one — on the r.h.s of the evolution equation, 
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namely an entirely new contribution to the JIMWLK Hamiltonian JfJJ). Obviously, these contribu- 
tions should not be considered a running-coupling correction, and we do not aim to compute them 
here. As explained in Sec. 13. II running-coupling corrections are usually identified by considering di- 
agrams with fermion-loop insertions, namely considering the formal large-TV/ limit with fixed Nfa s 
(the flavor expansion), and then making the substitution Nf — ► — 6/3o. However, in the present case 
this criterion is insufficient: similarly to running-coupling corrections, qq pair production appears 
at LO in the flavor expansion, as can be seen in Fig. [2J Moreover, the running-coupling correc- 





2tr(t a U Zl t b Ul 2 ) 
xt a U x ®t b Ul 



Figure 2: Virtual and real contributions at LO (a), (b) and NLO (c), (d) and (e). Vertically (as delineated 
by the background shading) they are grouped by the Wilson line structure induced by the target interaction. 
The latter are indicated on the diagrams via large solid dots, their analytical form is listed at the bottom. 
Diagram (e) induces a new structure via a qq loop that interacts with the target. 

tions and qq pair production are physically distinct only if the transverse separation of the pair is 
non- vanishing. Therefore, the identification of running-coupling corrections to the kernel requires a 
separation prescription. As wc shall sec below, the use of the dressed propagator in the background 
field automatically implements such a prescription. 

A related complication arises in the derivation of the JIMWLK Hamiltonian with running coupling 
because of the possibility of interaction between the background field and a produced qq pair. 
Dressing the virtual diagram, shown in Fig. [2 (a), poses no problems, since the gluon itself does not 
interact with the target 20 . For each LO virtual contribution in as well as their self-energy- 
like analogues, there is just one type of dressed virtual counterpart; e.g. in Fig. O diagram (c) 
corresponds to dressing of diagram (a). The eikonal interactions with the target field remain 

20 At this point one might be worried that some of intermediate quark-gluon vertices (whose location is to be 
integrated over) might lie on different sides of the target line - this is precluded by causality since the interaction is 
restricted to the light-like support of the target field 6+ = <5(x~ )j3{x). Even the dressed gluon line cannot propagate 
back and forth across this hyperplane. This will be borne out by the explicit calculation. 
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unaffected by the dressing. The dressed propagator in (c) reduces to the dressed propagator in the 
absence of the target field, and the running-coupling corrections may be directly obtained using 
the tools of Sec. 13.11 

On the other hand, real-emission diagrams, such as the diagram in Fig. [2(b), lead to two physically 
different generalizations upon inserting fermion loops; these are shown in diagrams (d) and (e), 
respectively. Diagram (d) contains the same eikonal factors as the LO diagram, while diagram (e) 
contains a new channel: the production of a well separated gg-pair, which interacts with the target 
field at two distinct points on the transverse plain z\ and 22. Instead of the single gluonic eikonal 
factor U z as in diagrams (b) and (d), in diagram (e) one encounters a factor 2tr(t a U Zl t b Ul 2 ). 
The distinction between the two, however, is lost if the transverse separation of the gg-pair (or 
the invariant mass of the pair) becomes small: in this limit the gg-pair in the final state becomes 
indistinguishable from a gluon. This (by necessity) is mirrored in the eikonal structure associated 
with the diagrams: 

lim 2tr(t a U zl t b Ul) = Uf . (39) 

Therefore, while in the virtual diagrams running-coupling corrections can be reconstructed solely 
from bubble sums as in Fig. [21 (c) , in the real emission diagrams these corrections are associated 
with both Fig. [21(d) and the local limit of (e). However, real and virtual contributions are related 
by the requirement of probability conservation: in the absence of interactions the dressed virtual 
terms must cancel the dressed real contributions exactly. In the following we show that it is 
possible to generalize the LO calculation of Sec. 12. II in a way that manifestly implements all these 
requirements. Specifically, the separation prescription between the running-coupling corrections 
that we compute and the qq production channel that we neglect, simply amounts to the replacement 
of 2tr(t a U Zl t b Ul 2 ) of diagram (e) by its local limit of l|39|l . namely U® b . As announced, this is 
automatically realized upon using the dressed propagator in the background field to compute the 
diagrams of Eq. (|10|) . 

The main ingredient in the LO calculation in [20] is the propagator of the fluctuation 8 A in the 
presence of the background field b + . The calculation sketched in Sec. 12.11 was performed in the 
A~ = axial gauge in order to reduce the number of diagrams that contribute. These diagrams are 
shown in Eq. I|10|) . The eikonal nature of the interaction in these diagrams implies that one only 
needs the "++" component of the fluctuation propagator, (SA+5A+) ab . This remains true for the 
dressed propagator. 

At LO, the gluon propagator can be expressed (see Eq. i|45J) below) in terms of an external ten- 
sor structure and the propag ator iG% b {x',y' ) of a masslcss scalar field in the adjoint color rep- 
resentation, propagating through the target field [3]. Our generalization will involve its massive 
counterpart iG^(x' ,y'). Since the structure of the propagator has been discussed extensively in 
the literature [3, 102], we only describe here the key ingredients that will be needed to apply the 
dispersive method and re-derive the JIMWLK equation with running coupling. 

Massive scalar propagator in the background field 

The starting point is a spectral representation of the massive scalar propagator 

Gt^v) (^ fc2 J l2 + tQ [<M^r, (40) 
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where we have combined a spectral integral over the virtuality J dk 2 , and the sum over all states 
at fixed virtuality, J d 2 k± dk~ / (2k + ), into one 4-dimensional integral J d 4 k. 

Here 4>1 c {x) are the solutions of the Klein-Gordon equation at virtuality k 2 = 2k + k~ — k 2 in the 
background field b + of Eq. J3J) : 

(-D[b} 2 -k 2 y b <f> b k c (x)=0 , (41) 

where a and b are color indices and c is a basis label. D^b] = <9 M — igb^ is the (adjoint) covariant 
derivative in the background field J3J. Explicitly, 

^ b{x)= / ^pe^S^-k 2 ) ) / rf 2 «- i ("-^[U;_ 1 i _ oo ( Z )] a6 , (42) 

reflecting the fact that the interaction with the background field can transfer transverse momentum 
from the target, while it cannot modify the conserved k~ component nor the vir tuality k 2 . The 
color structure in Ij42(l is carried by the adjoint eikonal factor [U~} _ (z)] ab defined by 



a u 



Pexp^ig / dz S(z )P(z) 
v~ 



(43) 



A straightforward calculation using l|42|) in i|4U|> yields: 



dk- 



2fc-(27r) 3 



[6{x- - y~)9{k-) - 6(y- - x~)e(-k-)] d 2 P± d 2 q ± 



-tp'X 



d 2 z± 



-i(p—q)z 



[u-} ty _(z)] ab e^y 



where the 4-momenta p and q obey the constraints 



(44) 



P 



2k- 



2k- 



Eq. I)44|l differs from the free massive propagator by the Wilson line Ij43(l that represents the inter- 
action with the target field. Since the interaction in l|43|l is localized at z~ = 0, this Wilson line 
reduces to U z (defined in J5J) if x~ > > y~ , to U\ if < < y- and to 1 otherwise. This 
implies [3, 102] free propagation 21 when the endpoints x^ and y- are both on the same side of the 
z- = hyperplane, while if they lie on opposite sides one encounters a three-step process: free 
propagation from the initial point onto the hyperplane, a current interaction with the background 
field, and free propagation again from there to the endpoint. Only in the latter case does the back- 
ground field induce a change of the transverse momentum in the propagator. It is important to 
note that these features carry over to the fluctuations propagator (8A+5A+) ab in both the massless 
and the massive case. 



21 With U = 1 there is no obstacle to performing the z integration in 1441 . 



21 



Massive gauge field propagator in the background field 



The relevant. "H — h" component of the fluctuation propagator at LO is expressed in terms of the 



scalar-field propagator iG^ b (x',y r ) as 



(6A+6A+) ab 



1 

IF 



D [B] x 4 Gq (x', y') + dl,iG a b (x',y')dl, 





' 1 " 


} 


9-. 



(45) 



Here an integration convention over 4-vector coordinates x' and y' is implied. Note that one may 
simplify the first term in this expression using the defining equation for Go, namely 

D 2 [b] x , Gf(x' >y ')=5^(x' -y>). 

We have chosen not to do this since Gq'{x' ,y') represents the propagating particle pole that will 
be affected by the dressing, while the remainder of the expression will not change. To illustrate 
that this is indeed the right procedure let us connect back to the free case, or alternatively restrict 
ourselves to virtual diagrams such as Fig. 0(c), by replacing 



D 2 [b] ab 



to obtain 



{5Ai5A+T b \ 



iGf(x,y) 



□ 



[x,y) 6 ab 



x'y' 



x'y' 



o- 



sab 



(46) 



(47) 



This coordinate-space expression may be readily Fourier-transformed and shown to coincide with 
the "++" component of the free axial-gauge propagator 22 of i|25b|) . which is (1 / k 2 ) 2k + / k~ . One 
can in fact relate the two terms in i|47[) respectively, to the contributions of longitudinal and trans- 
verse polarizations in the conventions of light-cone perturbation theory. To this end let us define 
transverse polarizations by 



A 



fc • 



-,e£,0) based on e T := (l,i(-l) A+1 ) /\/2 with A = 1,2. 



Here we used the notation (+,_!_,—)'' for the components of a four-vector. We may then recast 
Eq. (|25bp (before dressing the propagator), for our light-cone gauge (where n 2 = 0), as 



£2 I 9i* v 



n^k v + k^n u 
k ■ n 



(k ■ n)< 



fc2 T2 



: f_L A 



(48) 



\=1,2 



This uniquely identifies the second term in l|47() with transverse polarizations while the first term 
is related to longitudinal contributions. 

In P7|l the factors ^(x,y) correspond to the propagating particle pole 1/fc 2 in (|25 b>|> . which get 
multiplied by 1/(1 + n(fc 2 )) upon dressing the gluon. The point of separating the simple tensor 
structure of (1/k 2 ) 2k + /k~ by polarizations into two terms as in (|45|l and l|47(l is that the separation 



2 For our purposes the gauge vector n projects out the minus component of a vector: n ■ a = a 
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helps to extract the leading small-a; contribution: when used in the diagrams of Eq. , only the 
second term contributes a logarithm in 1/x, while the first term is sublcading. 

Eq. H47f) offers a straightforward route to apply the dispersive method for virtual diagrams such as 
Fig. |2](c). As explained above, in this case the dressed propagator does not interact with the target 
field, so the standard dispersive method described in Sec. 13. II directly applies. Following Eq. IpTTJl 
one replaces 1/k 2 by l/(fc 2 — to 2 ) or, in coordinate language 







I 




_n_ 


x'y' 


- □ + TO 2 - 


x'y 



After making this replacement in i|47|l the first term may be further split according to 



□ 



x'y' 



— W x'y' 



□ 



x'y' 



so the massive propagator finally takes the form 
1 



(SA+SA+^ ab 



V I m I free 



d- 



□ 



□ 



x'y' 



d- 



(49) 

5 ab . 
'(50) 

As we shall see below, the 5- function term does not generate a logarithm in 1/x, while the other 
two terms do, and will therefore be relevant for the JIMWLK equation. 

The dispersive representation of the dressing was derived in the free case and it therefore directly 
applies only to (|50|l and thus to the resummation of running-coupling corrections in the virtual 
diagrams in the JIMWLK Hamiltonian. Nevertheless, a similar procedure can be applied to real- 
emission diagrams, where it provides a natural definition of the running-coupling contribution, 
whose separation from the qq-p&ir production channel is a priori ambiguous. To this end we 
construct the massive propagator in the background field by reversing the replacements ()46f) in (I5U|) : 



d- 



l0 x'y' 



to 2 iG<£(x',y')+cP xl iG a J;(x',y')cP y , 





' 1 " 


} 


9~. 



(51) 



By using this propagator in both real and virtual diagrams, the real-virtual cancellation mechanism 
is in place, just as at LO. Thus, probability conservation is guaranteed. It is by this fundamental 
principal that the dressing of the virtual diagrams by fermion-loop insertions, which is uniquely 
determined by l|5U|l using the dispersive method, essentially dictates the structure of running- 
coupling corrections in the JIMWLK kernel as a whole. 

Generalization of the dispersive approach 

To incorporate running-coupling corrections in the derivation of the JIMWLK equation, we will 
modify the fluctuation propagator in analogy with Eq. 136f) . namely 



IT 



(SA+6A+ 



+ \ab 



1 

To 



dm 



T pv(m 2 ) (6A+6A+)% - (SAi6A+) a b 



(52) 
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where (<L4+<L4+)^' is given in J5TJ. Note this modification can be traced back to a substitution in 
the scalar propagator 

oc Snab , ,_a s f d 4 fc [Mx)4(y)] ab f d 4 k aY(~e~ l Q) [<M^l(y)] Qb 



G o (*.l0- n j (27f)4 k2 + l0 ~> 7(2^)4 8 tt fc2 + l0 (53) 



1 Z 1 dro 



A) 



2 



Pv (m 2 ) G a Jt(x,y)-G a b (x,y) 



where in the final expression we used the dispersive representation of the running coupling (|34|) . 
Restoring the tensor structure then yields (15211 . 

As explained above, our procedure, which uses the propagator (|51|l for both virtual and real- 
emission diagrams guarantees probability conservation. Its diagrammatic interpretation at the level 
of fermion loops can be read off Fig. [2 in the absence of interaction, (the limit U — ► 1), the sum of 
diagrams (d) and (e) reduces to diagram (c) and their cancellation is complete; as the interaction 
is turned on they become distinct and evolution takes place. At the same time a new channel, the 
production of a qq pair in diagram (e), opens up. Our procedure uses the leading- Nf corrections in 
the virtual diagrams to identify (define) running-coupling corrections. This amounts to a specific 
separation of the leading-A^ corrections in the real-emission diagrams into ones that are part of 
the running coupling and ones that constitute the new production channel of diagram (e). By 
virtue of Eq. ()39)l the local limit of diagram (e) is included in the running-coupling contribution. 
As a consequence, real-virtual cancellation holds separately for the running-coupling corrections, 
and the remainder, namely the new qq production channel can be computed separately: it is a 
well-defined NLO contribution to the kernel that is not associated with the running coupling. This 
is confirmed by an explicit calculation of the diagrams using light -cone perturbation theory [103]. 

Finally, let us briefly comment on the separate terms in Eq. 1)51)) . We begin by observing that, 
as in the LO case, the first term in 1)51)1 . — iS^ yl , is suppressed at large p~ and will therefore not 
contribute to the evolution equations. The two terms containing G m , however, will generate d~ 
contributions in the numerator, and will therefore contribute to logarithmically-enhanced terms 
at small x through J dp~ jp~ — ln(l/x). Next, we note that the two surviving terms are very 
different in nature: the last term is already present in the LO fixed-coupling calculation, while the 
middle term does not contribute there, as it vanishes in the m — > limit. It starts contributing at 
NLO. This is in keeping with the different origin of the terms: the transverse partial derivatives 
in the last term reflect the Weizsacker- Williams field structure of the LO emission kernel 1155)1 
that is entirely driven by transverse polarizations. The middle term, in contrast, arises entirely 
from longitudinal polarizations that are absent at leading order. As such it has an entirely new 
dependence on the transverse coordinates not present at LO. The corresponding diagrammatic 
calculation in [103] confirms this structure and the association of the terms with transverse and 
longitudinal polarizations; the latter enter there (starting at NLO) via instantaneous contributions. 



3.3 Application to JIMWLK 

The dispersive method introduces running-coupling corrections by performing a direct replacement 
of the gluon propagator in the leading-order calculation according to 1)52)1 . Therefore, the gener- 
alization of the JIMWLK equation to running coupling essentially constitutes of recalculating the 
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LO diagrams of I|1U|) with the massive propagator of (|51|l instead of that of (|45fl . In the following 
we first present an explicit calculation of one of the diagrams, and then generalize and derive the 
running-coupling corrections to the Hamiltonian H[U] in © . 

Let us first recapitulate the ingredients: in ljTU|) we consider the calculation of the leading logarithmic 
correction to the propagation of a fast moving dipole, a qq-p&ti. As explained in Sec. |3 this object 
can be represented as a product of two Wilson lines along the classical trajectories of the quarks. 
We take the trajectories to lie along the minus light-cone direction at x + = 0. The pair is then 
characterized by the transverse locations x and y of the quark and antiquark, respectively. 

The first diagram 23 in (|10|l contains all the information we need: the other two can be deduced 
from it. It describes the leading ln(l/x) contribution to the exchange of a gluon between the two 
quarks in the presence of Weizsacker- Williams background field. The result (see [20,23]) involves 
the fluctuation propagator and the Wilson lines that represent the external quark and antiquark: 




:g 2 I dx-dy~(5At6A+) a ; 

[e(-x-)U x t a + 9(x-)t a U x ) ® (9(y-)Ult b + 9(-y-)t b u£), 



which, after some algebra involving the fluctuation propagator l|45|l (sec below) may be recast as 

= - £?| ln(l/x) jd 2 z K xzy {[U]f t a U x ® t b Ul + [tf]f U x t a <g> Ult b 

-5 ab [U x t a ®t b Ul+t a U x ®Ult b \y (54) 

The four terms on the r.h.s are in one-to-one correspondence with the four diagrams in Eq. 
and K. is a purely two-dimensional scale-invariant kernel, given by 

(55) 

where we used the notation for "parent-" and "daughter-dipoles" introduced in (JSJ. The last 
formula shows explicitly that this leading order kernel changes sign at a circle whose diameter is 
defined by the "parent." To arrive at the final expression in l|54|) with the kernel of (|55f) we have 
performed the x~ and y~ integrals, using l|45|) with the explicit formula for GQ b (x,y), given by 14411 
with m — 0. These integrations give rise to factors 2p~ /p 2 and 2q~ /q 2 , respectively. The minus 
momentum components in the numerator compensate the explicit d~ in the denominator in (|45|) . 
while the additional J Q dk~ /k~ in Gg b (a;, y) leads to the factor ln(l/x) upon imposing cutoffs in x. 



Note that the first term in the curly brackets in (|45|l reduces to — iS x , yl and it does not contribute 



to the equation since it is suppressed at large p~ . 

Eqs. H54|) an d l|55|l summarize the result for one of the diagrams in i|10[l used in deriving the 
JIMWLK equation [20]. As announced we will now use this diagram to demonstrate how the 

23 This quantity was called Xxy in [20]. 
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JIMWLK equation can be promoted to running coupling by means of the dispersive method. To 
this end we use the substitution llo2l) obtaining 




dm' 



Pv(m 2 ) 




\ 



(to = 0) 



(56) 



/ 



Recall that pv(rn 2 ) is a universal object, the discontinuity of the coupling on the time-like axis, 
defined in . To compute running-coupling corrections using Ij56(l we need to repeat the calcula- 
tion of the diagram, but now with the massive propagator of (|51|1 instead of the massless one (|45|l . 
Doing this we find that Eq. (|54|) is unmodified while K XZ y is replaced by its massive counterpart: 



JC' 



d 2 pd 2 q c lpi - x -^ e^-y) 
(2tt) 2 p 2 + to 2 q 2 + m 2 



(57) 



where the two terms d 3 x d y and to 2 originate in the last and the middle terms in Eq. I|51|l . respectively. 
As in the LO case, the — iS x V y , term in Eq. I|51|l does not contribute in the small-x limit. 

Let us now compute the integrals in l|57|) explicitly. Denoting the "daughter dipoles" respective 
lengths by := |r^| (i = 1, 2) we obtain: 



K-xzy ~^xz V T\mK\{r\m) r2mK\{rim) — to Ko{rim)K {r2m) 



(58) 



where K. XZ y is the LO kernel of (|55|) and Kq(x) and K\{x) = —dKo(x)/dx are if-Bessel functions. 
These functions depend only on the lengths of the vectors not on the angles — a direct con- 
sequence of the angular integration in 157|) . Angular dependence appears here only through the 
derivative d x d y , giving rise to the factor 1C xzy in the first term in l|58[) . just as at LO. 

The next observation is that upon applying the dispersive method to the self-energy diagrams in 
Eq. I|10|) . in analogy with (|56|l . one recovers again the LO result with the kernel of i|57|) . Thus. 
all the diagrams share the same massive kernel IC x l zy of Eq. i|58|) . Since also the dispersive inte- 
gral of Eq. H56J) is common to all three diagrams, running-coupling corrections to the JIMWLK 
Hamiltonian H[U] of ® appear through the replacement: 



-fC 



xzy 



xzy 



0o 



dm' 



pv(rn 2 ) F{r\m 1 T2m) — 1 



(59) 



where we defined the "characteristic function" J r (rim,r2m), following the general discussion in 
Sec. 13.11 by extracting the LO kernel: 



xzy =-K-xzy J 7 {r 1 m 7 r 2 m) = JC xzy T {r 1 m 1 r 2 m) + T (r 1 m,r 2 m) 
J- T (rim,r2m) = r\m Ki{rim) r2mKi(r2fn), 

T h (r\m, T2m) 



(60) 



2 r 1 2 r 2 2 , . , , 
to K (rim)K {r2m) 



V\ ■ T2 



and used the fact that lim T{r\m, T2to) = 1. Here T T and T correspond to the first and second 



m— »0 



terms in l|58|l . which in turn originate in the d x d\\ and to 2 terms is (|57f) . respectively. As mentioned 
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in Sec. 13.21 fsee Eq. (|51|l there) these terms are associated with transverse (T) and longitudinal (L) 
gluons, respectively, which explains the notation in (|60|l . As we see, these two contributions to K.™ zy 
are rather different in nature: the first can be naturally written as a product of the scale-invariant 
kernel JC XZ y times a function of the lengths of the "daughter dipoles" , while the second brings in 
an entirely new structure that does not share the angular dependence of the LO kernel. 

We also introduce a notation for the "effective charge" 



fl(riA,r a A) := 



1 



dm- 



Pvim 2 ) J 7 {rim,r 2 rn) — 1 



(61) 



which plays the role of a scale-dependent a s /7r in all our calculations: the substitution of 

M xzy := R(nA,r 2 A) JC xzy (62) 



— K, 



xzy 



in © promotes all the fixed coupling diagrams entering the JIMWLK equation to running coupling. 
The running coupling corrected JIMWLK Hamiltonian then reads 



H[U] 



Z7T 



(63) 



All other ingredients in the JIMWLK equation (JSJ and the definition of correlator averages J2J 
remain unchanged. Eq. 162(1 . and with it (|63|1 . is our key result that fully implements running 
coupling in the JIMWLK equation. 



3.4 Borel representation of the resummed kernel 

To systematically define the perturbative sum and analyze power corrections in the JIMWLK kernel 
■M X zy, it is convenient to write a Borcl representation of the "effective charge" R, in the form of 
Eq. ®: 

1 f°° (n 2 \~ u 

Mxzy = ICxzy R(r 1 A,r 2 A,) = — J duT(u)l — \ JC xzy B(u, r x /i, r 2 /i), (64) 

where the Borel function is related to the characteristic function T (see Eq. ijHHl) by 



5,, sm iru 
ea 



dm ( m 



mf 



T(r\m, V2m) — 1 



(65) 



Owing to the different nature of the transverse and longitudinal contributions to T ', it is convenient 
to deal with them separately also on the level of the Borel function. To this end we write 



B(u,r 1 n,r 2 fi) = B T (u,rxfx,r 2 fx) + B L {u,r 1 ^,r 2 n), 



(66) 



corresponding to the two terms in (|60|l . Below we present explicit expressions for these functions 
in both momentum and coordinate space. Using (|57|l and l)65|l we obtain the following expressions 
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as Fourier transformations from momentum space: 

K xzv B-(u, W ^)= - e ««™ I ^e— e- 



dm 2 



o 



(2tt): 

p • q p . q 



m 

TT (-1) 



(p 2 + m 2 )(q 2 + ra 2 ) p 2 g 2 



2 \ ~~ W / 2 \ — u 



dvd, eip . ri e . iq . r2 p_s_ 

(2tt) 2 p^ q^ q 2 — p 2 



„ „t , \ 5,.sin7rw f d 2 pd 2 q 4 _,„.*.„ 

K J (27T) Z 



X 



dm ( m 



o 

d 2 p d 2 q 



(j? J (p 2 + m 2 )(q 2 + ra 2 ) 



(2tt) 2 p 2 -q 



m 2 J [l 2 

2 



(67b) 



These expression assume a perhaps more familiar form if we restrict ourselves to one-loop running, 
where, with T(u) = 1, expanding under the Borel integral amounts to replacing the powers of the 
momenta k = p,q, i.e. (k 2 / fi 2 ) " by the corresponding geometric series according to 



k 2 \ a s ((i 2 ) 



,2 



7T 1 + W) ln(fea/p2) 



(68) 



In fact, using (|64|l and 1)68(1 . the expressions in (|67|l can be directly mapped onto the expressions 
derived diagrammatically in [103]. 

We note that in the purely virtual case of Fig. [5] (c) the result simplifies in the expected manner: 
there the z integral may be performed and, in the absence of interaction with the target, it sets 
q = p. Then the sum of the integrands of transverse and longitudinal contributions l(67|l yields 



— u 



from K. B T horn K B L 

Via (|68(1 . this corresponds to (l/p 2 )a s (p 2 e~%)/ir, which is the expected contribution of a gluon 
of transverse momentum p 2 that does not interact with the background field. Note that there is 
an important cancellation at NLO (and beyond) between the separate contributions of transverse 
and longitudinal gluons. On the level of the dispersive integrals used to define B T and B L this 
cancellation is even more obvious: for KJ£ zy (c.f. Eq. it takes the form 



p q 



(p 2 + m 2 )(q 2 + m 2 ) p 2 
which has exactly the same interpretation. 



(70) 
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Direct Fourier integration of the perturbative sum obtained through I|tj8f) is not well defined owing to 
the Landau pole. In this respect these all-order expressions are largely symbolic. Correspondingly, 
the Borel integrals of both l|67a|) and (|67b(l converge for u — > oo only when the momenta are larger 
than the QCD scale A. Thus, when we exchange the order of integrations and perform the Borel 
integration after the Fourier transforms, the Borel integral will no longer be unambiguous for any 
value of the coordinates. As we shall see explicitly below, the coordinate-space Borel function 
exhibits poles along the integration path, at positive integer values of u. This reflects ambiguities 
in summing the series, which are indicative of power corrections (see the general discussion in 
Sec. . Knowing the location of the poles and the parametric dependence of their residues on the 
hard scales involved, one can estimate the effect of non-perturbative power corrections. It is the 
advantage of the Borel method to expose these power corrections in such a clear and concise way. 

With this in mind we now perform the Fourier transformation of the Borel functions. The all-order 
expressions lead to infinite sums that can be recast in terms of hypcrgcomctric functions. The 
calculation is done using integral representations of the Bcssel functions. The integrals over the 
dispersive mass intertwine the parameter integrations so that it becomes necessary to use Mellin- 
Barnes representations to decouple them. The Mellin-Barnes integrals are done last and lead to 
infinite sums of residues. Details of the rather involved algebra are given separately for B T and B L 
in Appendix 1X1 and IbI respectively. 

To express the results in a compact manner we define 

r> := max {r 1 ,r 2 } ; r< := min{r 1 , r 2 } ; £ ■= — • (71) 

r> 

We present two versions for each of the results for B T and B L , as a series in £ 2 and in 1 — £ 2 , 
respectively. This is important for numerical implementations near £ = and £ = 1, respectively, 
and will be used in Sec. 16. II below. In addition, the first form conveniently displays the singularity 
structure and the functional form of the residues, while the second form is convenient for performing 
an expansion in u, in order to extract the perturbative coefficients. 

The coordinate-space result for B T is: 



T i r i ' r 2 I 4e 3 \ sin(7r«) f 
>C XZ yB (u, rifj,, r 2 /x.) = 2 2 2 9 s "(1 - u)T(-u) 

_ k r(2-u + fc)r(i -u + k) 



i' 2 \ ' >^ 

oo 



T(l + k)T(2 + k) 
(tf (2 - u + k) + *(1 - u + k) - *(2 + k) - *(1 + k) + ln£ : 



fc=0 

X 



(72a) 



which can be recast in terms of a hypergeometric series 2 Fi using 15.3.11 of [104] 

ri • r 2 / 4e~* \ sin(7ru) u(l - u)T 2 (-u)T 2 (l - u) 
= ~if^[~>iJ^j 7r r(2 - 2u) 



(72b) 



x 2 Fi 



1 — u, —u 
2-2u 



i-e 
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A useful alternative expression to (|72b(l can be obtained using the identity: 



2f 1 



1 — u, —u 
2 - 2m 



i - 



2 r 1 



1 — u, 1 — u 
2-2u 



The result for B L is: 



>C XZ y B L (u, r X [i, T2fi) = - 



sin(7ru) 1 / 4 e 3 



E 



fe=0 



r(k + i-u) 



T(k + 1) 

x (ln(£ 2 ) - 2*(fc + 1) + 2*(fc + 1 - u)) 



(73) 



(74a) 



which again can be expressed in terms of a hypergeometric function, this time using 15.3.10 of [104], 



sin(7ru) 1 /4e 



r 4 (i- M ) 

T(2 - 2m) 



2f 1 



/ 1 — u, 1 — u 






\ 2~2u 



(74b) 



The result is fully symmetric in r\ and r2, in particular all coefficients in an expansion in powers 
of u (or a s ) have this property. 

Our expressions for B T and B L can be put to use in several ways. The first two pertain to obtaining 
information on the perturbation theory: 

• A generating function for the perturbative expansion in powers of a s (/x 2 ) (at a fixed scale). 
To this end one expands the results for the Borel function around u = and then integrates 
over u in (|64(l term by term. Assuming one-loop coupling (T(u) = 1) this amounts to the 
replacement: 

,2\\ "+ 1 



(75) 



Knowing the analytic structure of the Borel function, which has a leading rcnormalon sin- 
gularity at u = 1, it is clear that the expansion: (1) does not converge, and (2) it is not 
Borel summable; A definition of the sum in (|64|) is required. Starting from the expansion, a 
natural prescription is to truncate the sum at the minimal term. A more systematic regular- 
ization, which has been found useful in several applications in QCD, sec e.g. [52,54,101], is 
the Principal Value (P V) of the Borel integral in l|64() . 



In the next section (Sec. we will perform the perturbative expansion of the JIMWLK kernel, 
study numerically its apparent convergence in the first few orders, and compare it to the PV Borel 
sum. We also examine the approximation of the sum by scale setting. These issues are revisited in 
Sec. where we examine the effect of the newly computed higher-order corrections on the evolution 
in the case of the BK equation. 

The third application goes beyond perturbation theory; we can extract information on the power 
corrections from the integrand: in (|64|1 : 
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• The renormalon poles at positive integer u in B T and B L can be used to infer the parametric 
form and significance of non-perturbative power corrections that are expected to affect the 
evolution kernel. 



This is discussed in Sec. [3] for a single evolution step and re-visited in Sec. when solving the 
evolution equations numerically. 



4 Perturbation theory 

Having computed the Borel function B(u) entering Eq. l|64f) . we have essentially determined the 
expansion coefficients for i?(riA, r 2 A), and thus for the JIMWLK evolution kernel, to all orders 
in the large-/3o limit. To get the coefficients one expands the expressions for B T (u, fi/x, r 2 n) and 
B L (u,ri/j,,r 2 fj,) in Sec. 13.41 under the integral in Qfi4[). 



B T (u,r 1 n,r 2 fi) = 1 + ^&^(ri^,r 2 /i)w n ; B L {u 7 nfj,, r 2 /z) = ^ ^(n/z, r 2 /i)u n (76) 

n— 1 n — 1 

and integrates over the Borel variable term by term , to obtain: 



M xzy = R( ri A,r 2 A)!C xzy =}C xzy ^p- U + f^nl(^^^j b n {r lt t,r 2 n)J , (77) 

with 

b n (rifi, r 2 fj,) := b^(rin, r 2 fj,) + b^(r lf i, r 2 fi) • (78) 

The purpose of this section is to study the effect of these perturbative corrections to the JIMWLK 
evolution kernel. 

Explicit expressions for the expansion coefficients will be presented in Eqs. H8()|) and (|82(l below. 
Before looking into the details, let us briefly recall what one expects on general grounds following 
the discussion in Sec. 13.11 

• The all— order sum is renormalization scale invariant: Eq. i|tj4f> is renormalization-scale 
inavriant. This means that also (|77|l shares this property, if it is formally considered to all 
orders. However, the choice of the expansion parameter would make any finite order partial 
sum scale dependent. 

• The series in non-summble owing to infrared renormalons, which amount to 
power— suppressed ambiguities: Going beyond the level of a finite-order partial sum, 
in an attempt to compute i?(riA,r 2 A) to all orders, one finds infrared renormalon ambigui- 
ties. Let us explain how they arise here: since B{u) in (|64|) has a finite radius of convergence 
(u = 1), b n in (|77|l are characterized at high orders by geometrical progression with no sign 
oscillation. With the explicit n\ growth in JZZJ), it is obvious that the series would not con- 
verge, but instead reach a minimal term, and then start diverging. An optimal perturbative 



'Since we are using one-loop running coupling, this amounts to making the replacement of Eq. 1751 . 
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approximation may be defined by truncation at the minimal term. This, however, is incon- 
venient since the truncation order may (and in our case does) depend on all the scales in the 
problem (in our case r 2 , r\ and r 2 ). A more systematic definition of the sum can be made 
using the Borel sum (JHOJ). This integral does not exist as it stands, since the poles appear 
at real positive values of u, on the integration path. One can define the sum by shifting the 
contour above or below the real axis, or by choosing the Principal Value (PV) prescription. 
The differences between these choices are power-suppressed in the hard scales: they scale as 
integer powers of A 2 r 2 , A 2 r 2 and A 2 r| (up to logarithms, see Sec.|3J). The size of this ambi- 
guity is obviously controlled by the residues of B(u) in (|SU|) and can therefore be explicitly 
computed. It is also similar to the magnitude of the minimal term in the series. Since these 
ambiguities are expected to cancel against non-perturbative power corrections, they provide 
important information on such corrections to the kernel, that are otherwise unknown. We 
will return to discuss these corrections in detail in Sec. [5] 



Our default choice for defining the all-order sum of the series is the PV prescription. Having made 
this choice, we will examine the convergence of the perturbative expansion in the first few orders, 
and demonstrate (Fig. the dependence of this (apparent) convergence on the region of phase 
space as well as on the choice of the expansion parameter. 

Technically, the implementation of the PV regularization in the present case is complicated by the 
fact that Eqs. (|72|) and l|74|) contain terms with double poles 25 in u. We cope with this in the 
standard technique by isolating the double pole contributions and using integration by parts. The 
explicit expressions used in defining R are given in Appendix [EJ they are in turn based on the 
explicit expression for the first- and second-order poles extracted in Appendices [U] and [D] 

The numerical result for the PV Borel sum is shown in Fig. To appreciate the qualities of this 
prescription as means of defining the all-order sum, it is useful to compare it with increasing- 
order partial sums. The expansion of the Borel functions can be conveniently obtained starting 



1 — u , 1 — u 
2-2u 



with expressions <|74bf) for B L and l|72b[) with i|73|l for B T , which involve 2F1 

a special case of the type E function expanded in [107] (see Eq. (4.29) there). It is convenient to 
define the dimcnsionlcss variable 

4 e -§-2-y E 

0:=-^, (79) 

which appears as the argument of the logarithms in the coefficients. The expansion for B T takes 
the form 



xzy B (u, nn, r 2 fi) = - 



lr 2 - (rf + rf) 



rrr. 



1'2 



00 

E 



(80a) 



where 



6f (n/x, r U i) = ln(n) - , (80b) 



25 The current example is by no means unique in this regard. The same occurs also in the absence of large target 
fields, e.g. for the well-studied example of the Adler function D(Q 2 ) = Air 2 dU(Q 2 )/dQ 2 , see [91,105,106]. 
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bJir-Lfx, r 2 fj,) = i (bj(nfj,, r 2 /x)) - y - ^— -j- 



> . 2 i(e 2 in(a) 2 1 



63 (ri/i,r 2/ u) = 



1-e 



(1-e 2 ) 2 i-e 

-1(1 -a ln 3 (0)-iln 2 (0)C 2 ln^ 2 ) 



:Li 2 (l-^ 2 ) 



(80c) 



£2 ,,2 



(c 2 + i) (m-f) in(a + Li 2 (e 2 ))-^ 



ln(fi) 



+ 2£ 2 ln 2 (e 2 ) ln(l-£ 2 )+3£ 2 ln(£ 2 ) Li 2 (£ 2 ) + (2 £ 2 + 2) Li 3 (l - £ 2 ) 
-2£ 2 Li 3 (£ 2 )- 



4 , 10^ C3 _ 4e2(1 _ a 



2 = l-£ 2 



dS 2 , 2 (z) / dz in (|80d|) is the first occurrence of a Nielsen polylogarithm, 



S a ,b( z ) 



(-1) 



(a - 1)!6! y £ 



^ ,_o-l 



In- 1 ^ ln°(l-£z) 



(80d) 



(81) 



At higher orders in the expansion one encounters [107] higher Nielsen polylogarithms, as well as 
other harmonic polylogarithms [108]. 

The longitudinal part, B L starts at order u, corresponding to O(f3oa 2 ). The expansion takes the 
form: 



n=l 

where, from (|74b|l . we find 

1 ln(£ 2 ) 



(82a) 



IC xzy &f(ri/x,r 2 /x) = - 



!C xzy b 2 (rin,r 2 fi) 



1 ln(£ 2 )ln(ft)+2Li 2 (l-£ 2 ) 



IC xzy &3 (r!/z,r 2 /x) 



' > 

1 1 



rl 1 - £ 2 



in(^ 2 ) n 2 (e 2 ) - 2 n 3 (e 2 ) + 2 c 3 - 4Li 3 (i - e 2 ) 



- ln(£ 2 ) ln 2 (fi) + 2 ln(fi) Li 2 (l - £ 2 ) 



(82b) 
(82c) 

(82d) 



It is interesting to note that there is a qualitative difference between the leading behavior of the 
transverse and the longitudinal contributions, respectively. The sign of the LO, transverse contri- 
bution to the kernel M xzy in 177|) can be directly read of Eq. 1)550: it is 



M > inside the circle: r\ + r 2 < r 2 



Mi 



* XZ y 



< outside the circle: rf + ri > r 



while the NLO longitudinal contribution to Ai xzy , in (|82b)l . is always positive. 



(83) 
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With the explicit expressions of (|80[) and (|82[1 we can now study the convergence of the perturbative 
series in (|77|) and compare the finite-order results to the PV definition of the all-order sum. Of 
course, at this point we must 26 make a choice of scale [i. Since M xzy = R{r\A, r 2 A) K, xzy is 
renormalization invariant (i.e. //-independent) one may choose any function of n, r 2 and r as long 
as the convergence of the perturbative series is good enough. As we shall see, this choice makes a 
significant difference for the apparent convergence at the first few orders. 

A priori, a natural choice (see Sec. 12.2(1 may be the "parent dipole" size: /i 2 = c/r 2 . A posteriori, 
knowing the coefficients in (|80J) and Q82[l. one can optimize the scale to reduce the size of subleading 
corrections. One possibility is to introduce BLM scale setting [80] by requiring that the entire NLO 
coefficient &i(ri/i,r2^) would vanish identically 27 : 

fcxzybiiriH, r 2 \i) + JC xzy bf(r 1 ^, r 2/ u) = 

2 4 f 5 e ln(£ 2 ) / 2r 2 \ I 

=► mLm = ^ exp | -- - 2 JE + S-J^l (l+ r2 _ r2 >_ ^ j j ■ (84) 

In our case, however, this approach is too simplistic: since the leading-order kernel vanishes on the 
circle r\+r 2 = r 2 , as implied by l]83fl . while the NLO contribution associated with B L does not, 18411 
develops an artificial singularity within the parturbative region. Consequently, with this particular 
choice of scale, higher-order terms would not be bounded on the circle, while no special difficulty 
would arise there otherwise. While better possibilities for optimizing the scale exist — for example, 
setting the BLM scale such that the NLO transverse contribution would vanish bj(ri/i, r^ix) = 
- we will not take this avenue here. Let us just note in passing that the multi-scale nature of 
the problem, which is reflected in the dependence on the parent and daughter dipoles — or in 
momentum space by the separate Borel powers of the transverse momentum before and after the 
interaction with the target (Eq. I|67|) ) — renders any scale-setting procedure in terms of a single 
coupling unnatural. Further progress in this direction will be reported in Rcf. [103]. Here we 
consider a couple of simple examples for the scale: 

4 e -5/3-27 B g e -5/3-27 E 

it = ^ and u, 2 = ^ , — . (85) 

r 2 r 2 + rf v ' 

The choice of the constant 4e~ 5 / 3 ~ 27E f» 0.24 is motivated by Eq. lfM|t . 

In Fig. |2| we show a numerical comparison of the order-by-order expansion with the PV Borel sum. 
As expected, at short distance scales perturbation theory is well-behaved 28 : the first few orders 
in (|64J) gradually approach the PV Borel sum. For larger distances, the corrections become large, 
and the minimal term in the series is reached sooner (see Fig. 01(b)). Eventually, for ?'iA of order 
of a few times 10 _1 , the series diverges right from the start. Note that even in that region the PV 
Borel sum is uniquely defined, but as we shall see in the next section, power corrections are large 
as well. 



26 To avoid the ultraviolet problems of a conventional fixed-coupling calculation one in fact needs to choose fi as 
some function of T\ , T2 and r instead of a single homogeneous scale independent of the configuration encountered [45] . 

27 To express the YZ X zy in the transverse part in terms of dipole lengths only we substitute: —r\ ■ T2 = (r 2 — r\ — 
r|)/2, as in 1551 above. 

28 If we allow the constant in 1851 to vary significantly compared the choice indicated, the convergence in both 
Figs. |3] (a) and (b) becomes appreciably worse, an obvious consequence of the logarithmic terms in the coefficients, 
sec JSSJ and 1551 . 
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Figure 3: R(r\A, r 2 A) as defined with a PV regularization of the Borel sum 



and the convergence of 



perturbation theory at the first few orders, according to ¥171 . with fi 2 — 4 exp ( — | — 2-je) /r 2 (in (a)) and 
(i 2 = 8 exp ( — | — 27b) / (rf + r|) (in (b)). R is shown as a function of r\, with rA = 1CP 4 , r || ri \\ r 2 and 
r 2 = r 1 + r. 



It is interesting to further observe that there is a significant difference in the apparent convergence 
in the first few orders between Figs. (a) and (b): the "parent dipole" scale setting in Fig. 0(a) 
has significant corrections even at short distances, while with the scale choice of Fig. (b) the first 
few orders provide a better approximation. It should be emphasized though that Fig. shows a 
particular situation where r 1; r 2 3> r. If the dipole sizes are of comparable magnitude, there is no 
significant difference between these two choices of scale 29 . The fact that "parent dipole" running 
fails to approximate the kernel for r±, r 2 3> r, makes a difference for evolution when R s 3> n, r% S> r. 
We observe that parent dipole running in this region generically underestimates the perturbative 
sum. Accordingly, in Sec. 16.21 we will see that evolution is slower with the "parent dipole" ansatz 
than with the Borel sum. 

In Fig. 0] we compare the PV Borel sum with two ad hoc running-coupling formulae, the "parent 
dipole" running and the "square-root daughter-dipolc" running: 



1 



1/2 



parent = < S (c 2 /^) /tT = ^ ln (c 2 /(r 2 A 2)) ? R ^ = K ( C >?) «f (c 2 /r 2 ) J A- 

(86) 

In these expressions we fix the scale c 2 = 4e~3~ 2 7-E i as motivated by Eq. 1)85(1. This ensures that 
these couplings would match the PV Borel sum soon enough when both r 2 , r\ <C 1/A. This is shown 
on the right panel of Fig.0J We note here that c 2 here is much smaller than the value c 2 — 4 arising 

29 For clarity, let us comment that we have chosen to show the function R with r\ parallel to r2. This avoids the 
apparent divergence (r i ■ ) _ 1 in the definition of R L ; this divergence is exactly canceled in the BK kernel when R 
is multiplied by K xz , y . 
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Figure 4: Comparison of ad hoc implementations for effective couplings R calculated with the principal 
value (PV) prescription, the "square root" prescription [a s (c 2 /r 2 ) a s (c 2 /r 2 )] 1//2 /n, and the parent dipole 
running a B (c 2 / V 2 ) / V . On the left panel the size of the parent dipole is fixed, rA = 0.1 (and r-2 = ri +r); on 
the right panel the ratio of the dipole sizes is fixed, r\ = 2r, r^ = Sr. 

definitions of running coupling can be made to match reasonably well at small enough dipole sizes, 
at least in parts of phase space, if the scales are chosen appropriately. On the other hand, at small 
dipole sizes, already at 0.1-0.3 x 1/A the coupling diverge strongly. This related of course to the 
smallness of the effective scale at which the coupling is evaluated, /j^g ~ 4e _ 3 _27E /r 2 . As we shall 
see below, the effect of this divergence on the evolution rate (in the BK approximation) is somewhat 
less important than for R. This is to be expected, as evolution rate is an inclusive observable, and 
as such has reduced sensitivity to the phase-space region where the effective coupling is large. 



5 Non-perturbative power corrections 



After discussing the perturbative features in some detail, let us now turn to estimates the non- 
perturbative corrections to the JIMWLK kernel M. xzy . The renormalon poles in R allow us to 
determine the type of power corrections, and their parametric form, based on the singularities of 
B T and B L . The original expressions in l|72|) and (|74|l have simple and double poles at integer values 
of u. As explained in the Appendices O and iDl we perform integration by parts over the double 
pole terms to arrive at a representation of the Borcl integral where the integrand is composed of 
simple poles only; in this way the double poles have been converted into a simple ones accompanied 
by a logarithmic-enhancement factor, ln(fi) where O is given in Eq. (|79|l . In this representation, 
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the explicit expressions for the renormalon poles are the following: 



IC xzy B T (u,r 1 fi,r 2 ^) =fC 




T(l -m + n) 



u — m T(n + 1) I r 2 /i 2 J T(n — 2(m — l))r(m — 1) 



EE 



(1 _£2)" 2(-l) m ~ r ' 

r(n +1) (u — m) 



dm 



4e~3 



r(2m- (n + 1)) 



r 2 u 2 / r(?7i)r(m — l)r(m — ri)r(?7!. + 1 — n) 



regular contributions 



(87a) 



and 



OC OC 



EE 

m— 1 n—m 



-(-l) m 1 /4e"3\ (1-C 2 )™ (T(n + l-m)) 



oo m — 1 

EE 

m— 1 n— 



2(-l)™-™ )™ 1 ^ 

u — m r(n + 1) r> dm 



r(n + 1) (r(m)) 2 r(n + 2 - 2m) 
T(2m- (n + 1)) 



4e"3 \ 

r 2 u 2 J (r(m)) 2 (r(m-n)) 2 



regular contributions 



(87b) 



To determine the ambiguity associated with each renormalon we should first sum transverse and 
longitudinal contributions, and then isolate the residue at fixed u = m. 

Power corrections are expected to follow this ambiguity structure. Let us therefore introduce a 
non-perturvative parameter C m , of order 1, relating each power correction to the corresponding 
renormalon residue. The underlying assumption is that genuine non-perturbative effects would 
be of the same order as the ambiguities. The parametric dependence of the corrections on the 
hard scales then follow directly from the residues in Eq. I|87|) : they are written as powers of r 2 A 2 
times some function of £ 2 , with additional logarithmic terms in fl. The dependence soft scales is 
subsumed in the coefficients C m . In doing so we write the full kernel M. xzy as a sum of perturbative 
contributions M^ zy and power corrections SAixzy'- 



M xzy = M p J y 



OO 



SM {m) 

/ j J l xzy 

rn — 1 



As an example we explicitly present the power correction SA4xl y corresponding to the leading 
u = 1 renormalon: 



In (£ 2 ) - 2 In 



4e~f ~ 2lE 
r 2 A 2 



- r? A 2 ea 



(89) 



37 




Figure 5: Rpy and the expected range of power corrections for fixed rA = 0.1 (left panel), and fixed ratio 
ri = 2r, r-2 = 3r (right panel) with ri parallel to ri- The central solid line corresponds to the principal 
value result. The shaded regions estimate the relevance of power corrections by adding and subtracting 
n\residue\(u = m). The inner band takes the first residue at u = 1, and the outer hashed band is the sum 
of absolute values of contributions from all residues. The dashed line shows the result when the u-integral 
has been cut at M max = 0.75, before the first pole is encountered. 



Fig-0shows R and the expected size of power corrections as functions of n, for fixed parent dipolc 
r = 0.1 /A, and for fixed ratio of dipolc sizes r2 = r±/2 (here ri = ri+r and r2 is parallel to r±). The 
power corrections are assumed here to be of the order of the ambiguity in choosing an integration 
contour in the Borel plane, ir times the absolute value of the residues at u = m = 1, 2, 3 . . .. The 
relative importance of renormalon poles can be seen from the width of the error band in Fig. |SJ 
The contribution from the pole at u = 1. Eq. (|89|l . is shown separately. It is evident that this first 
pole strongly dominates, except very close to 1/A, as expected from the hierarchy of powers in the 
analytical expressions. 

The key feature of the power corrections in Fig. [S] is the fact that they quickly die away at small 
scales 30 . This feature is very robust and it does not depend on the detailed prescriptions used 
to estimate the size of the power corrections. Just to give an example, one might estimate the 
uncertainty by cutting off the u integral at some value it max < 1, before one encounters the first 
renormalon pole. This is shown in Fig.0for u max = 0.75 with dashed lines. Clearly this prescription 
leads to an alternative estimate for the pcrturbative sum for R which is roughly within the error 
band of the full PV integral. 

From Fig. El it becomes strikingly clear that the poles limit the precision of our knowledge of the 
(non-perturbative) evolution kernel at large distances. We will see in Sec. El when we discuss 
evolution in the context of the BK equation, that for Q s near 1 — 2GcV, where large dipolcs 
contributes significantly to the evolution, power corrections must be taken into account if one 
wants to quantitatively determine the evolution. If we wish to apply JIMWLK or BK evolution 

30 Note that on the left panel of Fig.|S]r2 approaches r = 0.1/A as n — > 0. 
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starting from Q s -values in this range, both the evolution rate and the generic functional form of the 
initial condition for the dipole cross section N Yo ^ xy receive sizeable non-perturbative contributions. 



6 Evolution in the BK approximation 



Our discussion so far focused on the kernel of the JIMWLK equation. We computed and resummed 
running-coupling corrections to the kernel and extracted an estimate for non-pcrturbative power 
corrections. All these affect only the "effective charge" in front of the leading-order kernel, although 
in a way that strongly depends on the size of the evolving dipole, and that varies significantly over 
the transverse phase space of the newly produced dipoles. In order to explore the consequences 
of these results on the evolution, one should clearly consider the solution of the equation over a 
sufficiently large range in rapidity. Only in this way would it be possible to examine the sensitivity 
of observable quantities, such as the rate of evolution, to the corrections computed. 

The purpose of the present section is to perform a first exploration of this kind, by considering 
the evolution in the case of the BK equation. To this end, let us repeat the derivation of the BK 
equation, as described in Sec. 12.11 starting with the JIMWLK equation with the Hamiltonian of 
Eq. which includes running-coupling corrections. A simple way to see the way the corrections 
enter is to first write the leading-order BK equation with the kernel separated according Eq. 1101 
into exchange and self-energy-like diagrams: 



d Y N YiXy = ^ / 



aJj^]_ r usiy 2 ) ^ Usin 2 ) ^ 

^ f^xzy i^xzx i^yzy W^-V 

7T 7T 7T 



X 



(Ny,xz + N Y , Z y - Ny.xy - N Y , XZ N Y ,. 



The result obtained using l|t)3|l by repeating the steps of Sec. 12.11 simply amounts the substitu- 
tion (|62[) . separately for each of the three terms in 190|) : 

—it xzv -> M xzy = 2 R{rxK, r 2 A) K, xzy - i?(nA, n A) K, xzx - R(r 2 A, r 2 A) JC yzy (91) 

where the vectors t\ arc defined in (jSJ). Our final result for the resummed BK equation is then 
written as 

d Y N Y , xv = ^ Jd 2 z M xzy (n y , xz + N Y , zy - N Yxy - N Y , XZ N Yzy J . (92) 



6.1 Numerical implementation 

The numerical solution of the BK evolution equation l|92|l requires care. In the following we briefly 
describe the choices we made in our implementation. Firstly, we consider only translationally- 
invariant and spherically-symmetric solutions; we set x — > and N xy — > N\ y \. Despite this, the 
^-integral in Eq. (|92|l renders the evolution equation in essence two-dimensional. The simplest 
way to perform this integral is to discretize the two-dimensional space using a finite regular square 
lattice. This was used in [45] in order to compare BK with JIMWLK evolution. This approach, 
however, restricts the ratio of ultraviolet and infrared cutoffs (the ratio of the lattice size to the 
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lattice spacing) to, at most, ~ 10 4 , and hence strongly limits the Y range over which the shrinking 
correlation length R S (Y) can be resolved on such a lattice. A solution to this problem is to use 
discretization with higher resolution power at small distances, as was done in [77]. 

To achieve this we discretize N r on an even logarithmic scale, r n = r exp(riA), using ~ 250 points 
with r m i n = e~ 22 /A and r max = 1/A. The two-dimensional integral in Eq. l|9"2|) is evaluated using 
nested Simpson integrations in the (log|z|, argz) coordinates. While the values of \z\ and \y\ are 
restricted to discrete revalues, \z — y\ is not. Thus, it is necessary to interpolate the value of 
N\ z _ y \ in Eq. I|92(l using the known points N Tn . We have checked the stability of the simulation by 
changing the number of discretization points, and verified that this leads to a negligible change in 
the results. 

In order to evaluate the R- functions in the kernel M. of 1)91(1 . we need to integrate over the Borcl 
function B(u). However, as described in Sec.^ one cannot directly use the expressions in Eqs. (|72ajl 
to (I74b(l when integrating over the positive real axis, because of the presence of single and double 
infrared renormalon poles there. As described in App. [EJ the double poles were converted analyti- 
cally into simple ones by means of integration by parts. The integration over the simple poles was 
performed by a numerical implementation of the Principal Value prescription. This implementation 
is based on different analytic formulae, depending on the value of £ = r < /r > : for £ 2 < 0.8 we use 
<|123a[> and (|123b|) for R T and R L , switching over to ifT^jl and lfT23|) when £ 2 > 0.8. The integrals 
are completely stable when the threshold value of £ 2 is changed within the interval (0.01, 0.99), and 
our default value of 0.8 is chosen to optimize the evaluation speed. The numerical evaluation of 
the kernel is further stabilized by pulling the integrands of the three i?-functions appearing in 1)91)1 
under a single integral. 

If the parent dipole is much smaller than the daughter dipoles, r <C r\ ~ T2, the kernel M. contains 
large terms that cancel each other almost completely, leading to numerical difficulties. This can 
be dealt with by recognizing that the large contributions arise in fact from the leading terms in 
the series expressions of Borel functions Eqs. (|124fl and 1)125)) appearing in 1)92)1 (in this parameter 
region 1 — £ 2 RiO). Combining the first terms analytically stabilizes the evaluation of the series. 

The evolution in rapidity Y is performed with the second-order Runge-Kutta method, with a "time" 
step of 5Y = 0.05 . . .0.1. Again, we observed no significant difference in the results upon varying 
the step size. 

6.2 Simulation results 

The key feature of JIMWLK and BK evolution is scaling of the solutions with Q S (Y) after the 
initial-state effects have died out. This scaling is exact in the fixed-coupling limit and has been 
shown to be retained to a very good approximation in simulations with ad hoc implementations of 
running coupling. Also our all-order resummation result shows this feature, as is demonstrated in 
Fig. El 

The initial state in our simulations is N- U1 i(r) = 1 — exp[— r 2 /rp], where rg sets the initial scale. The 
evolution tends to flatten N(Y, r) until the scaling shape has been reached. However, the shape 
remains considerably steeper than with fixed coupling, see Fig. E|(b). From this plot we also see 
that the scaling shape is quite similar for the running coupling corrections computed in this paper 
and the ad hoc prescriptions shown. Thus, the asymptotic shape of N(Y, r) is largely insensitive to 
the functional dependence of the scale of the coupling on the phase space. 
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Figure 6: (a): The evolution of the function N(r) as a function of Y = Inl/x, shown in intervals of 
AY — 4. After initial settling down the function evolves towards smaller r while approximately preserving 
its shape. The main effect of the running coupling is to slow down the evolution at small r. (b): The shape 
ofN(Y,r) using the computed running -coupling corrections (PV) as well as some running-coupling models, 
taken at Y where the saturation scale R 3 (Y) — 0.005/A, where we define R 3 through N(Y,R 3 ) — 1/2. 
Different running-coupling forms display similar shapes for N(r), albeit slightly steeper for the parent- 

_ 2 / 2 

dipole running; much shallower than the initial state N- m \ = 1 — e r /r °, but much steeper than the fixed 
coupling shape shown for comparison. 



In Fig. [7| (a) we show the evolution of the saturation scale R S A as a function Y; here the full PV 
Borel sum result is compared with the parent-dipole running and the "sqrt" ansatz of l|8(j[l . The 
correlation length R s characterizes the scale where N(Y, r) is changing and saturation sets in; here 
we define R S {Y) through the condition N(Y,R S (Y)) = 1/2. The initial R S (Y = Y ) is w 0.1/A. 
We see that R S (Y) rapidly evolves towards smaller values. 

We note that a different definition of R S (Y) would naturally give somewhat different curves. How- 
ever, when the system has evolved sufficiently long so that it reaches the scaling form, all reasonable 
definitions should yield R S (Y) values that differ only by a constant factor. Above all, the evolution 
speed 

1 dQ 2 s (Y) 1 dR 2 s (Y) 

A( Y) ~ Ql{Y) dY - ~RW) dY ' Qs ~ 1/Rs (93) 

should be the same. 

As with the ad hoc prescriptions, the most interesting consequence of running coupling on JIMWLK 
and BK evolution is the fact that it drastically slows down the evolution compared to the purely 
fixed-coupling case by restricting the active phase space to within one order of magnitude of the 
characteristic scale Q S (Y) [45]. This qualitative feature does not depend on the details of the imple- 
mentation of the running-coupling effects. The precise value of the evolution speed, as expressed by 
the evolution rate X(Y) on the other hand, does depend on these details. For fixed coupling, scaling 
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Figure 7: (a): the evolution of the saturation scale Rs(Y) with the computed running-coupling correc- 
tions (PV), in comparison with a few different running-coupling models. The evolution starts at initial 
R 3 (Y — Yq) w 0.1/A. The horizontal dotted line shows the value of R„ where N(r) was shown in Fig. [(J 
(b): the evolution speed X(Y) plotted against the saturation scale R 3 (Y). The evolution starts at large R a 
values; the rapid increase at R a A ~ 0.1 is an initial-state effect. Several choices of scale are presented for 
the "parent dipole", demonstrating the sensitivity to the scale and the fact that upon tuning the scale, this 
model too can be made to agree with the resummed result for X(Y) at small R s . 



with Q s is exact, and A becomes a constant proportional to a s (/x 2 ) [79]. With running coupling, 
scaling is approximate, and A becomes a function of Y that will receive both perturbative and 
non-perturbative corrections via i?(riA, ^A), which we already examined numerically in Fig.JSJ 

While using Eq. I|93|l with N(Y, R S (Y)) = C (where C — 1/2 was used above) is certainly possible, 
during the initial stages of the evolution it is sensitive to the particular value of the constant C 
chosen. A more robust definition is the one used in [45], using a 1/r 2 moment of the evolution 
equation as an operational definition of A: 

N f d 2 r f 

A(F) = ^ / — / d 2 z M xzv (N(Y, r\) + N(Y, r\) - N(Y, r 2 ) - N(Y, r\)N(Y, r 2 2 )) (94) 
We will follow this convention here. 

Our results for \{Y) ar e shown in Fig. [7| (b). Here we plot A against R S (Y); in this way the 
result is independent of the initial value for Y . The initial-state effects arc visible as a very rapid 
increase of A near the initial (large) R s . As soon as these effects die out, the evolution approaches 
an initial-state-independent curve (see Fig.[8J). On this curve N(r) has reached the scaling shape 
shown in Fig. El and A decreases slowly as R s decreases at large Y. This is entirely driven by the 
logarithmic decrease of the coupling at short distances. 

In this plot one can also observe the sensitivity of the evolution rate to running-coupling corrections. 
We show the resummed result computed by the PV prescription (full line) in comparison with LO 
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evolution where the scale is set in variety of ways. LO parent-dipole running with fj? = A/r 2 
strongly under-estimates the evolution rate even at very small correlation lengths R s . This clearly 
demonstrates the significance of higher-order running-coupling corrections, which are included in 
the PV Borel sum. 

Interestingly, as seen in Fig. [3(b), one can approximate the PV Borel sum result (at sufficiently 
small R s ) with a variety of different functional forms, provided that an appropriate choice of scale 
is made. Thus, despite the large differences in the "effective charge" between the resummed result 
and the models considered — which are particularly significant for large dipoles — the differences 
in A, at sufficiently small R s , reduce to an overall phase-space independent multiplicative factor 
in the scale of the coupling. With the specific scales given in Eq. I|85[l the "square root" ansatz 
approximates the PV Borel sum well, whereas parent-dipole running slightly underestimates it. In 
the latter case, choosing /i 2 = 2e3 _27E /r 2 yields a very good approximation to A for R S A< 10 -2 , 
but somewhat over-estimates it at larger R s . The possibility to approximate A in l|94|) well as 
a function of the correlation length (provided it is sufficiently small) in terms of a phase-space 
independent coupling, such as the parent-dipole one, is probably related to the scaling property of 
N(Y, r), namely the fact that its shape, which determines the weight given to different final states 
in (|94|l is invariant. 




Figure 8: The evolution speed X(R S ) for different initial conditions Ni n i(r) = 1 — exp(— r 2 /rg), using the 
PV Borel sum. As the initial-state effects vanish, X(R S ) approaches a universal evolution trajectory. 

In Fig. [8] we show the evolution rate A(V) as a function of R s using different initial conditions. The 
trajectories approach a universal curve, which is independent of the initial condition. 

In Fig. (a) we show the evolution rate as a function of R s using the perturbatively expanded 
kernel, Eq. (|77|) at different truncation orders. Here the scale of the coupling is chosen to be 
/i 2 = 8e~ 5 / 3 ~ 2TB /(r 2 + r|) as in H85|) . We note that in contrast with the Borel sum, the perturbative 
expansion suffers from a Landau pole when the daughter dipoles become too large. These Landau 
singularities obstruct the phase-space integration, so the large region must be cut out. To avoid 
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Figure 9: (a): Evolution using perturbatively expanded kernel, Eq. I|77|l with fj 2 = 8e _5 ^ 3_27E /{r\ + r|). 
(b): Estimating the non-perturbative influence on X(R S ) by adding and subtracting n x (residue at u — 1) 
to the PV evolution kernel. The effect of the residue dies out quickly as R s becomes smaller. 

any visible effect by this cut we performed the simulation of Fig. 0(a) with relatively small initial 31 
R s . In the small R s regime shown, one observe good convergence towards the PV Borel sum. At 
sufficiently small R s only the leading-order result differs significantly from the PV one. This is of 
course specific to the (almost optimal) choice of scale made here. Much larger corrections appear 
for a generic scale choice as we have already seen in Fig. [7| 

Fig. El (b) shows the effect of the estimated power corrections on evolution. The relative importance 
of the power corrections decrease with Q S (Y): the active phase follows Q S (Y) towards the ultravi- 
olet. With increasing Q S (Y) the uncertainty induced by power corrections dies away very rapidly. 
Despite this general trend, it is apparent that a quantitative determination of evolution speed for 
current experiments where Q s = 1 — 2 GeV, requires to take power corrections into account. To 
test the predictive power of the approach, it is important to compare several distinct observables, 
and address the universality of the relevant power corrections. 

7 Conclusions 

We have presented a first calculation of running-coupling corrections to both JIMWLK and BK 
equations. We have shown that both these equations, which were originally derived at leading 
logarithmic accuracy with strictly fixed coupling, can indeed be promoted to the running-coupling 
case: the general structure of the equations remains the same, while the kernel receives corrections. 
Running-coupling corrections are singled out from other radiative corrections in the following way: 

31 To study the fixed— order perturbative result at larger R s with this or a similar phase-space dependent scale one 
would have to consider the effect of the cut or introduce another rcgularization of the Landau pole. 
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with or without running coupling, the r.h.s of the evolution equation of any correlation involves 
just one more Wilson line (the produced gluon) than the l.h.s (the evolving object); other radiative 
corrections entering at the NLO involve up to two more Wilson lines. The number of additional 
Wilson lines grows further at higher orders. It remains for future work to fully generalize the non- 
linear small-a; evolution equations to NLO. Some work in this direction, which is complementary to 
ours, has already been done [56,57,76]. A full NLO generalization exists of course for the linear case 
of the BFKL equation, and our expectation is that it exists in the non-linear case as well. Since the 
BFKL equation is known in full at NLO, one can make a useful comparison. This goes beyond the 
scope of the present paper. One important fact, however, that we do learn from the BFKL limit, is 
that running-coupling corrections constitute a significant part of the total NLO correction. This, 
we expect holds also in the JIMWLK and BK cases. 

The significance of the running coupling in the context of JIMWLK and BK evolution has been 
acknowledged long ago. For one thing, in fixed-coupling evolution the active phase space extends 
way into the ultraviolet, while practically any dependence of the coupling on the scales present at 
a single evolution step reduces the active phase space to within one order of magnitude around 
the correlation length R s [45]. For this reason all simulations of JIMWLK or BK evolution have 
been performed with running coupling, implemented using some educated guess as to its scale 
dependence. The simulations presented in this paper are the first where the scale dependence 
is computed from QCD perturbation theory. While qualitative features of the effects of running 
coupling on the evolution — e.g. the decrease of the evolution rate X(Y) with decreasing correlation 
length R s (y) at large Y — are similar to what has been observed before, the scale factor itself 
has been determined here for the first time. As shown in Fig. [7| using the correct scale is crucial 
in obtaining a quantitative estimate of the evolution rate. With this predictive power, we can 
hope that non-linear evolution equations with running coupling would become directly applicable 
to small-a; phenomenology at LHC energies. 

Technically, running coupling corrections are computed, as usual, by focusing on a specific set of 
diagrams where a single gluon is dressed by fermion-loop insertions, making use of the linearity of 
Po in the number of flavors Nf. To perform this calculation we utilized in this paper the dispersive 
approach, where the all-order sum of vacuum-polarization insertions is traded for a dispersive 
integral in the "gluon mass". While this technique has been used extensively in the past, using 
it to compute the diagrams entering the JIMWLK Hamiltonian is non-trivial for two reasons: 
first, the large- Nf limit by itself is not sufficient to disentangle running coupling effects from 
the new production channel of a qq pair; second, the presence of the strong Weizsacker- Williams 
background field may interfere with the "dressing" by interacting with the fermions. Guided by 
the correspondence between real and virtual corrections owing to conservation of probability, we 
could nevertheless disentangle running coupling corrections from other contributions and generalize 
the dispersive technique to the case of the background field. To this end we derived a dispersive 
representation of the dressed gluon propagator in the background field. This formal development 
may well have other applications. 

By computing the JIMWLK Hamiltonian using the dispersive technique we could go beyond 
the NLO O(j3oal) running coupling corrections to the kernel, and resum C(/3q _1 q;™) corrections 
to all orders in perturbation theory. Besides the obvious advantage of obtaining an exactly 
rcnormalization-scale invariant kernel, the all-order calculation offers a unique window into the 
non-perturbative side of the problem. The fact that infrared-finite evolution equations can be 
established in the high-energy limit, does not mean of course that the dynamics governing the 
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evolution is purely perturbative; non-pertubative dynamics affects the evolution though power- 
suppressed corrections. By performing an all-order resummation one can get access to this infrared 
sensitivity by looking at the ambiguity in separating perturbative and non-perturbative contri- 
butions; these are the infrared renormalons. Using Borel summation we identified explicitly the 
ambiguities in defining the perturbative sum, and in this way established an estimate for the para- 
metric dependence the of non-perturbative effects on the hard scales involved and the potential 
magnitude of these effects. 

We find that both perturbative and non-perturbative corrections modify the evolution kernel in 
a non-trivial way. In particular, these corrections depend on all the different scales present: the 
"parent dipole" r as well as the two "daughter dipolcs" , n and r-i . In this way different final states 
are weighted differently. An interesting feature at NLO and beyond is the appearance of two classes 
of contributions: one that is proportional to the LO kernel, which is associated uniquely with trans- 
versely polarized gluons, and one that is not, which is associated with longitudinal polarizations. 
Both perturbative and non-perturbative corrections to the kernel, propagate through the evolution 
and affect the rate at which the saturation scale Q S (Y) flows towards the ultraviolet with increasing 
energies. Conversely, the saturation scale itself determines the active phase space, leading to the 
decoupling of soft modes at large Y. 

We studied here the effect of the newly computed corrections on two levels: first by looking at 
the "effective charge" controlling the contribution to the r.h.s. of the evolution equation, and 
then by solving the BK equation numerically and studying the effects on the evolution itself. We 
found, on both levels, that the running-coupling corrections are significant. Our simulations have 
shown that non-perturbative corrections are strongly suppressed at high energies, reflecting the 
fact that the increase of Q S (Y) with energy moves the active phase space along with it towards 
shorter distances. At presently accessible energies, where the saturation scale is estimated to be 
Q S (Y) ^1 — 2 GeV, power corrections arc definitely relevant. Going from high to low energies, 
one may view the breakdown of the perturbative evolution as the onset of the Soft Pomeron. If 
we accept this premise, our calculation opens a new way to think about the Soft Pomeron: as the 
energy is lowered the number of power-suppressed terms that need to be included in the evolution 
kernel increases. A possible determination of the power terms from data can thus be attempted 
at intermediate energies, before the power expansion breaks down. A detailed examination of this 
idea must address the distinction between the initial condition to the evolution and the corrections 
to the kernel itself. 

In conclusion, we made an impotent step in extending the framework of non-linear evolution equa- 
tions at small x to include running-coupling and power corrections. Nevertheless, this endeavor is 
by no means complete: small-a; evolution at high densities presents many challenges, some of which 
we touched upon in this paper. This includes for example the generalization of the non-linear evo- 
lution equations to the NLO; the detailed comparison with the NLO BFKL kernel; understanding 
the relation with DGLAP evolution and higher-twist corrections to the twist expansion, gaining 
better understanding of the initial condition for the evolution; the development of (experimentally 
accessible!) obscrvablcs of different degree of inclusivencss, which are sensitive to the dynam- 
ics underlying the evolution; and the determination to non-perturbative corrections affecting the 
small-a; evolution. 
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A B in coordinate space 

We begin with T T in d dimensional coordinate space. The regulator is only needed to have sep- 
arately finite terms in (|97|l and may be removed once the cancellation of infinite contributions is 
manifest by setting d — > 2. The simplest way to obtain a dimcnsionally regulated expression for 
T T is to start from Eq. I|57|) and to promote the momentum integrals from 2 to d dimensions. One 
then exponentiates the denominators using two Schwingcr parameters t\ and i 2 - This leads to 

d OO OO 

r l r lY J dt r (« e -- 2 -4^ J ds s -m + D e -s^-^pi (95) 





^ ?T (r 1 m, r 2 m) 



and a d dimensional generalization of Eq. (|60J) 
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From here we obtain B T according to (|65|) : 
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We will first perform the m 2 integral by using ^El 
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The T — 1 structure serves to regulate at m = 0. To capture its effect it is expedient to make use 

-m 2 (*i+t ) *l+*2 2 

of =5 — — = — J da e~ m a to decouple the m 2 and ij integrals 

o 

d OO OO 



r 2 (|) V 4 4 



o 

*l+*2 ^ 

/ da [°° dm 2 [ '-^ ) 



» 



(/i 2 )"a"- 1 r(l-«) 



The m 2 integral converges only if Re(u) < 1. u = 1 will be the location of the first renormalon 
pole. At this point we may set d = 2 since all expressions are now explicitly finite: 

= -(»r r j r jjdh e-^Jdt 2 it e-* fr+^rq-tQ . (99) 



One might guess that the poles of T(l — u) reflect the location of the renormalon poles, but this is 
not true: once one restores the pre-factor to recover B, T(l — u) is actually canceled by sin(7rw)): 

5 / o\ U — 00 00 



R T/ \ es» M 2 ^r 2 r 2 V f (d / 2 +i) -M/fl / , . ,- W 2+i) 

£r (u, rim, r 2 m) = — — - — — / dU ti e «i I dtot y ' 

y ' ' ; r(u + i) U 4^ J 1 7 2 



e 



x(ti+t a r (100) 

which would be completely regular if not the t integrations expose poles at positive u. 

This can be uncovered by using a Mcllin-Barnes representation for the {t\ + i 2 ) u factor in this 
expression (see [109] for a textbook on applications of the Mellin-Barncs technique in field theory). 
We write 

_1_ r r(-u + u)T(-u) t\ 
2m J r(-u) t 2 



(*x + w = ^ / d » v z , (ioi) 



7 

where the path 7 connects —77 — zoo to 77 + zoo (77 real and infinitesimal) in such a way that it 
separates all left from all right poles of the integrand. [To recall the definition: poles originating 
from factors T(a + v) and T(a — v) arc called left and right poles respectively The left hand side is 
recovered as a series in either t\/t% or i 2 /ii by summing residues. The sign of ln(ti/t2) determines 
which way to close the contour and in turn the form of the series by selecting the relevant residues.] 
This decouples the t integrations which provide additional v dependent T functions allows us to 
write B T as 

— 77 — \~ ti OO 

(r 2 /4) M 1 f /r 2 \' y r 2 (l-u + 7/)r 2 (l 



B 1 {u,r ll i,r 2l x) = -e^ - — / dv -| — > > '- . (102) 

1 (u + 1)1 {— u) 2m J \ r 2 J —u + v v 



— r] — ZQQ 



The integral can now be done by summing residues. We have to distinguish two cases: the sign 
of ln(? ,2 /r2) determines where to close the contour. It proves useful to emphasize the symmetry of 
this procedure by expressing the result via r< , r> and £ 2 = J' 2 - / ?' 2 : 
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ln(rf/7"2) < 0, closing to the right: 
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ln(rf/r|) > 0, closing to the left: 
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(n + l)(n + 1 - u) 

This is identical to (fTU3|) . 

Using shift identities for the digamma functions, 2-0(71 + 1) = i/j(n + 1) + tp(n + 2) — ^-j- and 
2^(71 + 2 — m) = ^>(n + 1 — u) + i/'(?t. + 2 — u) + ra+ }_ M , we arrive at our final result 
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(105) 



B B L in coordinate space 

The calculation of the coordinate expression for B L is slightly simpler than that of B T and uses 
the same techniques. 

We start from (|65|) for the longitudinal contribution and note that its contribution to the m = 
subtraction term vanishes: 

r nf , . 5,shi7ra [°° dm 2 ( ' m 2 \~ u _ z . . 

K xzy B L {u, ri ^r 2 ^)^~^ u / — — K%%. (106) 

n Jo TO V M / 
Abbreviating /(u) = _ e f u 2i^ anc i 

using the parameter the appropriate parameter representation 



for the K factors this turns into 
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For the remainder of the caclulation we present only minimal comments, 
-separate off the m integral- 



oo oo 







-the m integral requires Re(it) < 1- 



(2tt) 2 

-use MB to factor (ti + t 2 ) 
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-separate the t integrals and perform them (Re(— d/2 + v) < — 1 and Re(— d/2 + u — 1) — v < —1) 
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-rearrange to expose convergence of the MB integral- 
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-evaluate MB independently for both cases r 2 > r\ and r\ > r 2 , find sums of residues to agree, set 
<i = 2 to obtain the final result - 
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C Pole and renormalon structure of B T 

To expose the pole structure of B T we find it easiest to start from l|72b(l and rewrite it as 



B T (u,r 1 fJ l ,r 2 fJ.) = [ i^-r \ - yT/TTT-i — 2F 



r \ " sin(7r U ) < l ~ u ) ( r (! - U ) T (- U )) 
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f Mt/4e-lV" 1 Tjl-u + n) 

— gr^Tiy^J i^iw^R r(2 - u)r( - u+n) (m) 

Here all is regular but the factor T(2 — u)T(n — u) which will cause at most double poles at positive 
integers u = m=l,2,.... In the following we will encounter Laurent expansions around given pole 
locations and use the abbreviation e = u — m. 

• Double poles arise via combined divergence in T(2 — (m + e))T(n — (m + e)). For a fixed pole 
location u = m = 2, 3, . . . , 00 the second factor appears to contribute as long as n < m. The 
Laurent expansion around u = m in these cases can generically be written as 
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We note that coefficient of the double pole also vanishes for n = m so that the sum of the 
double pole terms may be expressed as 
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• As is obvious from the Laurent expansion in i|113|) we inherit single poles at the double pole 
locations. Additional single poles arise when 

— only the first of the factors contributes a pole, i.e. for n > m > 2, (the coefficients here 
vanish, however, where n < 2(m — 1)) 

— or where only the second of the factors does contribute, i.e. for 2 > m > n (this is 
limited to m = 1, n = 0, 1): 



51 



The Laurent expansion around any of these additional single poles is given by 

Adding all these leads to the following sum of single pole contributions which mark the power 
correction to to B T 
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dm I \ r 2 /x 2 / r(m)r(m — 1)T(to — n)T(m + 1 — n) 
As a power correction each of the terms at fixed m comes with an a priori unknown coefficient. 



D Pole and renormalon structure of B 1 
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Renormalon poles again appear at u = m = 1, 2, 3, . . ., with everything regular but the factor 

r(i-u)r(i -u + n). 

• Double poles at m arise where both L(l — (m + e))T(n + 1 — (m + e)) diverge, i.e. for 
1 < m > n + 1 . The Laurent expansion around these pole locations it = m is 
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so that the double pole part of the longitudinal contribution may be written as 



y double pole 



EE 
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(119) 



As for B T , we inherit single poles at the double pole locations. Additional single poles arise 
when only the first factor in T(l — (m + e))T(n + 1 — (to + e)) diverges, i.e. for 1 < to < n + 1. 
The Laurent expansion around these new pole locations reads 
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Some of the residues vanish (where 2to > n + 2). 

Combining the two contributions leads to an expression for the single pole part of B L of the 
form 
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E Principal value definition of the perturbative sum 

The small £ 2 series differ strongly from the 1 — £ 2 series so that we will deal with them separately. 

• Closer inspection of (|64|) via (|72a|l and l|74b|) reveals, that the double poles in the obstructing 
terms have their origin in a u derivative of a T function, i.e. they arise from differentiating 
simple poles. To make this manifest we rewrite (|72a|l and (|74at as 
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The double pole contributions arise from factor — u + n) and as such they can be 

converted into single poles under the Borel integral by a partial integration in the expression 
for R, Eq. Ij64(l . For this, we note that boundary terms of the partial integration in (jf>4|> vanish 
to display the results separately for R T and R L as 
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(123a) 



and 
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• The single pole expressions given above converge only very slowly near £ = 1. To circumvent 
this we also present single pole expressions based on l|72b(l and (|74b(l . For the transverse 
contributions we write 
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- double pole 

where B T is taken from (|72b(l and the double pole part of the integrand is derived in Ap- 
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pcndix[U] It is given by 
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The first term can now be integrated with a principal value prescription (it only contains 
single poles) while (for one loop running) the last term is given a meaning via integration by 
parts and gives only a boundary contribution. To this end we recall that for one loop running 
we may set T(u) -> 1, T'{u) -> and identify In (y 2 / 'A 2 ) = g^;n?n : 
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double pole 



oo m—X 

EE 

m— 2 n—2 



4e"3 

TCa 2 



(l-f) n 2(-l)' 



r(2m- (n + 1)) 



r(n + 1) m r(m)r(m — l)T(m — n)T(m + 1 — n) 



We have checked that this result matches with a very costly evaluation of R T via (|123f) . R L 
can be treated analogously 
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K 



xzy 



^2 ) B L {u,r 1 n 1 r 2 y) 



oo m — 1 

EE 



- double pole m=l n=0 

r(2m- (n+1)) 



2(_i)"»-» (1 -£ 2 )™ 1 / 4e"t 
(tt-m) 2 r(n+l)r§. I r?,A 2 



(r(m)) 2 (r(m-n)) 2 



(125b) 



55 



and to one loop accuracy 
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The expressions on the r.h.s. of all the Borel integrals shown in this section are in general well 
defined with any arbitrary contour connecting to 00 in the complex plane that avoids the poles 
of the integrand on the real axis. All such expressions share the same perturbative expansion (be 
they partially integrated or not!). The uncertainty introduced by the this freedom of choice in the 
contour is qualitatively determined by the size of residues of these poles. An important consistency 
condition for using the partially integrated version to study non-perturbative corrections is therefore 
that the residues of the integrand remain the same as in the original version l|64[l . This is a generic 
consequence of the fact that T(2 — u + fc) (the function contributing the poles) has no logarithmic 
terms in its Laurent expansions. Appendices iGl and iDl contain an explicit discussion of the pole 
part of B T and B L . 
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